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ABSTRACT 


The  search  for  a  universal  solution  of  the  equations  of  motion  for  a  satellite  or¬ 
biting  an  oblate  planet  is  a  subject  that  has  merited  great  interest  because  of  its 
theoretical  implications  and  practical  applications.  The  discovery  of  such  a  solu¬ 
tion  should  motivate  a  reassessment  of  both  the  theories  that  exhibit  singularities 
and  the  physical  effects  implied  by  singularities.  The  practical  importance  of  such 
a  solution  is  the  efficiency  of  simple  analytic  formulas  in  predicting  simultaneously 
the  paths  of  large  numbers  of  satellites  in  a  multitude  of  orbits.  Here,  a  complete 
first  order  solution  to  the  problem  of  a  satellite,  perturbed  only  by  the  oblateness 
of  the  Earth,  is  displayed.  The  orbit  is  free  of  singularities  for  all  parameters  and 
is  valid  for  1000  revolutions  with  a  relative  error  of  the  order  J 2  10-6. 
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a  semi  major  axis  of  the  initial  instantaneous  ellipse 

C  COS  *0 

e  eccentric1,1)-  of  the  initial  instantaneous  ellipse  (0  <  e  <<  l) 

h  cons’  mt  value  of  ^  along  an  orbit,  approximately  equal  to  the 

initial  magnitute  of  angular  momentum 
i  inclination  of  reference  plane 
t0  initial  value  of  i 

J2  oblateness  coefficient  of  the  planet  (coefficient  of  the  second  harmonic  in 
the  expansion  of  the  gravitational  potential) 


p  h2 )GM  where  G  is  the  gravitational  constant  and  M  is  the  mass  of  the 
planet 

R  equatorial  radius  of  the  planet 

r  radial  distance  from  the  center  of  the  planet  to  the  satellite 
s  sin  z’o 
t  time 

1 0  initial  value  of  t 
v  * 

r 

\  gravitational  potential 

<5  latitude  of  the  satellite 

<p  longitude  of  the  satellite 

6  angle  from  line  NON'  to  the  satellite  measured  in  the  reference  plane 
(Fig.  5.1) 

6i  initial  value  of  6 

ft  longitude  of  line  NON'  (Fig.  5.1) 

ft  initial  value  of  ft 

^  argument  of  perigee  of  the  initial  instantaneous  ellipse 
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I.  INTRODUCTION 


A  characteristic  feature  of  practical  satellite  orbit  prediction  is  that  the  engi¬ 
neer  may  deal  with  numerous  satellites  in  a  great  variety  of  alternative  orbits.  Un¬ 
der  these  and  many  other  such  circumstances  analytic  relations  which  can  quickly 
approximate  an  orbit  may  be  far  superior  to  large  numerical  models.  While  many 
analytic  methods  have  been  developed  for  the  artificial  satellite  age,  most  are  not 
used  in  practical  orbit  prediction  because  they  violate  one  or  more  of  the  following 
principles: 

•  The  method  should  provide  a  solution  that  is  significantly  more  accurate 
than  the  two-body  solution. 

•  The  real  physical  effects  of  the  orbit  should  be  easily  distinguishable  in  the 
solution. 

•  The  solution  should  be  universal;  it  should  be  valid  for  all  orbital  parameters. 

The  motivation  for  this  research  was  the  desire  to  develop  a  method  for  satellite 
prediction  that  would  embody  these  characteristics. 

In  this  analysis,  a  solution  to  the  equations  of  motion  of  a  satellite  around  an 
oblate  planet  is  found  by  use  of  a  variation  of  the  perturbation  technique  known 
as  the  Method  of  Strained  Coordinates.  The  orbit  is  valid  for  1000  revolutions 
with  a  relative  error  of  10“6.  The  solution  which  is  valid  for  all  eccentricities  and 
for  all  inclinations,  was  obtained  by  extensive  use  of  the  symbolic  manipulation 
program  MACSYMA. 

The  analysis  begins  with  a  background  discussion  of  some  of  the  competing 
satellite  orbit  theories.  There  is  then  a  development  of  the  equations  of  motion 
beginning  with  a  derivation  of  the  two-body  solution.  The  various  forces  which  act 
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to  disturb  the  two-body  orbit  are  highlighted;  a  more  thorough  discussion  is  given 
for  the  effects  of  oblateness.  There  is  a  complete  treatment  of  the  perturbation 
technique  as  the  equations  of  motion  are  solved  in  detail.  The  complete  first  order 
solution  is  displayed  as  a  function  of  coordinates  and  as  a  function  of  the  orbital 
elements.  In  addition,  an  independent  analysis  of  the  polar  and  equatorial  orbits 
is  performed  to  serve  as  a  check  of  the  general  solution. 
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II.  BACKGROUND 


A.  INTRODUCTION 

The  theory  of  flight  of  artificial  satellites  is  closely  related  to  classical  celes¬ 
tial  mechanics,  one  of  the  oldest  and  most  highly  developed  branches  of  science. 
The  equations  which  describe  the  motion  of  an  artificial  satellite  axe  in  principle 
identical  to  the  equations  of  motion  of  natural  celestial  bodies.  It  is  not  surprising 
then  that  the  results  originally  derived  in  classical  celestial  mechanics  have  been 
freely  used  to  explain  the  motion  of  artificial  celestial  bodies. 

The  foundations  of  classical  celestial  mechanics  were  established  in  the  eigh¬ 
teenth  century  when  Clairaut,  d’Alembert,  Laplace,  Lagrange,  and  Euler  intro¬ 
duced  theories  and  analytical  methods  to  explain  the  large  deviation  of  the  Moon 
from  an  elliptic  orbit  due  to  solar  attraction.  These  theories  all  supplemented 
or  complemented  the  pioneer  work  that  had  been  done  by  Newton.  Newton  had 
correctly  indicated  the  Moon’s  variation  in  eccentricity  and  inclination  and  the  re¬ 
gression  of  the  nodes;  however,  his  published  theory  accounted  for  only  half  of  the 
motion  of  perigee.  (A  century  later  an  unpublished  work  was  found  to  contain  the 
full  explanation.)  Clairaut  in  1749  was  on  the  verge  of  substituting  a  new  law  of 
gravitation  for  the  Newtonian  law  when  he  found  that  second  order  perturbations 
removed  the  discrepancy  in  the  motion  of  perigee. 

In  that  same  century,  Euler  investigated  and  developed  the  perturbative 
function  and  began  the  development  of  the  method  of  variation  of  parameters, 
which  was  later  extended  by  Lagrange. 
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Lagrange,  Laplace,  and  Poisson  all  advanced  the  discipline  through  their 
investigation  of  the  stability  of  the  solar  system.  Later,  in  the  nineteenth  cen¬ 
tury,  the  use  of  Hamiltonian  mechanics  was  used  to  great  advantage  by  Delaunay 
whose  work  influenced  the  innovators  of  artificial  satellite  orbit  prediction,  namely 
Brouwer,  Kozai,  and  Garfinkel. 

The  introduction  of  the  modern  computer  and  the  launch  of  the  first  artifi¬ 
cial  satellites  had  a  profound  effect  on  the  science  of  celestial  mechanics.  Classical 
celestial  mechanics  had  been  essentially  a  contemplative  science  with  the  principle 
aim  being  to  study  the  laws  of  motion  of  existing  heavenly  bodies.  In  contrast,  the 
science  of  the  flight  of  artificial  satellites  is  an  active  engineering  science  concerned 
with  determining  or  predicting  relatively  short  term  orbits,  and  in  many  cases  con¬ 
trolling  the  satellite’s  motion  through  on-board  propulsive  devices.  It  was  in  the 
exciting  climate,  following  the  successful  launchings  of  the  first  artificial  satellites, 
that  most  of  the  new  methods  of  satellite  orbit  prediction  were  developed. 

B.  THE  USE  OF  HAMILTONIAN  MECHANICS  IN  SATELLITE 

ORBIT  PREDICTION 

There  are  certain  schools  of  thought  in  dynamic  astronomy  and  theoretical 
physics  that  support  the  loyal  use  of  Hamiltonian  mechanics  [Ref.  1,  p.  228], 
and  many  of  the  new  methods  in  general  perturbations  take  advantage  of  the 
elegant  formalism  offered  by  the  Hamiltonian  method.  The  Hamiltonian  method 
is  referred  to  here  as  the  formal  process  of  writing  the  equations  of  motion  for  a 
satellite  in  the  canonical  form: 


dqr  _  dH 
dt  dpr 


dpr  _  -dH 

dt  dqT 


(r  =  1,2...  3n) 


(2.1) 


where  qr 


generalized  coordinates 


( H  is  the  Hamiltonian  and  L  is  the  Lagrangian,  L  =  T  —  V,  kinetic  energy  - 
potential  energy). 

The  solution  to  (2.1)  may  be  written  down  if  a  function  S  can  be  found,  where  S 
is  any  complete  solution  of  the  Hamilton-Jacobi  equation 


It  should  be  noted  that  (2.2)  is  tractable  only  if  the  variables  qr  and  t  are  separable 
within  5  and  H . 

Although  this  present  analysis  does  not  use  Hamiltonian  mechanics  to  solve 
the  equations  of  motion,  a  summary  of  its  use  is  merited  since  the  advances  in  this 
area  have  been  an  invaluable  contribution  to  general  perturbation  theory. 

In  celestial  mechanics,  one  may  formulate  a  Hamiltonian  that  represents  the 
gravitational  attraction  of  the  central  force  and  add  to  it  terms  which  are  “per¬ 
turbing  Hamiltonians.”  Various  sets  of  canonic  variables  may  be  chosen  with  the 
goal  of  expressing  a  zero-order  Hamiltonian  in  a  simple  form  and  the  higher  order 
effects  in  an  iterative  fashion.  Each  term  then  may  be  dealt  with  through  a  succes¬ 
sion  of  canonic  transformations.  Delauney  introduced  a  systematic  procedure  for 
isolating  parts  of  the  Hamiltonian  and  then  generating  a  suitable  transformation 
in  successive  steps.  A  particular  feature  of  his  approach  is  that  a  periodic  term  in 
the  Hamiltonian  may  be  eliminated  with  each  canonical  transformation  [Ref.  2]. 
While  DeLauney  used  a  procedure  that  eliminated  one  term  at  a  time,  von  Ziepei 
devised  a  technique  that  eliminates  one  angular  variable  with  each  transformation. 
This  method  reduces  the  numtxr  of  degrees  of  freedom  while  at  the  same  time 
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imparting  to  the  transformed  Hamiltonian  a  symmetry  of  the  unperturbed  system 
[Ref.  1].  The  eliminated  variables  are  referred  to  as  ignorable  coordinates  since 
they  do  not  participate  in  the  solution  of  the  transformed  equations  of  motion  but 
can  be  recovered  after  the  solution  has  been  obtained. 

Using  von  Ziepel’s  method,  Brouwer  devised  one  of  the  most  notable  general 
perturbation  theories  [Ref.  3].  Prior  to  Brouwer’s  method,  all  previous  work  in 
artificial  satellite  theory  had  written  the  Hamiltonian  as  a  Fourier  series  in  the 
mean  anomaly  with  coefficients  that  were  infinite  series  in  powers  of  the  eccen¬ 
tricity.  Brouwer  used  an  elliptic  approximation  for  the  potential  and  obtained  a 
complete  first  order  theory  with  some  second  order  development  using  canonical 
transformations.  The  essence  of  Brouwer’s  method  was  to  write  equations  (2.1)  as 


dLi_dH_  dfy  _  dH 

dt  dl,  dt  dLj 

where 


Li  =  VGMa  h-M 
L2  =  Li's/1  -  e 2  l2  —  u> 
L3  =  L2  cos  i  u  =  n 

are  the  Delauney  variables  where 


a  -  semimajor  axis 

i  -  inclination 

u)  -  argument  of  perigee 


e  -  eccentricity 
M  -  mean  anomaly 
fl  -  longitude  of 

ascending  node 


(2.3) 


Referring  back  to  the  Hamilton  Jacobi  equation  (2.2),  the  function  S  is  used  as 
a  generating  function  to  find  a  new  Hamiltonian  that  leads  to  simplified  canonic 
equations.  By  choosing  5  correctly,  Brouwer  was  able  to  find  a  canonical  transfor¬ 
mation  from  the  Delauney  variables  to  a  set  of  double  primed  Delauney  variables 
(L'\  /”)  such  that  (2.3)  have  the  form 
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dL" 

~r-=0 

dt 


dt 


dH’ 

dLj 


(2.4) 


where  H *  is  the  Hamiltonian  expressed  in  terms  of  double  primed  orbital  elements 

by 


L\  =  \/GMa"  l”  =  M" 

Lj  =  L"\/l  -  e"2  /"  = 

Lj  =  L![  cos  t"  /J  =  n" 

The  double  primed  orbital  elements  are  related  to  the  unprimed  elements  by 


a  =  a"  +  ea  e  =  e"  +  ee 
i  =  i"  +  ei  M  =  M"  +  eM 

w  =  w"  +  ew  n  =  n"  +  f  n 

where  the  quantities  ea,ee,ei,eM,tw,eCl,  are  periodic  functions  of  the  double 
primed  orbital  elements. 

Equations  (2.4)  are  solved  for  the  double-primed  Delauney  variables  which 
can  then  be  expressed  in  terms  of  the  original  variables. 

The  results  initially  obtained  by  Brouwer  were  not  valid  at  the  critical  incli¬ 
nation  of  63°. 4,  and  they  were  questionable  when  either  the  inclination  or  eccen¬ 
tricity  were  near  zero  [Ref.  4].  O.  K.  Smith  devised  a  method  for  dealing  with  the 
problem  of  zero  inclination  and  eccentricity  [Ref.  5],  but  Brouwer  subsequently 
challenged  the  validity  of  his  method.  Later,  Lyddane  [Ref.  6]  was  able  to  suc¬ 
cessfully  remove  the  restrictions  on  small  values  of  eccentricity  and  inclination  by 
reformulating  Brouwer’s  work  in  terms  of  an  alternate  set  of  variables. 

Brouwer  used  the  central  force  term  as  the  first  approximation  for  the  po¬ 
tential  since  there  is  no  exact  solution  to  the  equations  of  motion  of  a  satellite 
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under  the  influence  of  the  more  complex  potential  described  by  an  oblate  planet 
having  axial  symmetry.  Other  authors  have  attempted  to  introduce  a  new  poten¬ 
tial  which  approximates  the  Earth’s  potential  better  than  the  central  force  term 
alone  and  also  leads  to  an  exact  solution.  Most  notable  has  been  the  work  of 
Sterne  [Ref.  7]  and  Garfinkel  [Ref.  8].  Sterne’s  potential  function  accounts  for 
most  of  the  standard  potential  through  the  second  harmonic  and  it  leads  to  a 
solution  of  canonical  constants  that  are  free  of  first  order  secular  perturbations. 
The  remaining  effects  of  the  earth’s  oblateness  and  other  forces  are  allowed  for 
in  the  perturbing  Hamiltonian  which  causes  the  cix  canonical  constants  of  the 
unperturbed  solution  to  undergo  variations  with  time.  Sterne  provided  the  inspi¬ 
ration  for  Garfinkel’s  method  which  is  essentially  the  same  as  Sterne’s  but  more 
developed.  Garfinkel  included  the  second  and  fourth  harmonics  and  arrived  at 
a  solution  that  is  reducible  to  quadratures.  Garfinkel’s  original  solution  was  not 
valid  at  the  critical  inclination,  but  in  a  later  paper  he  removed  the  singularity 
through  a  variation  on  his  perturbation  theory  [Ref.  9]. 

While  Garfinkel  did  his  analysis  in  spherical  coordinates,  Vinti  [Ref.  lOj, 
derived  a  potential  expressible  in  oblate  spheroidal  coordinates.  In  his  original 
analysis,  Vinti  introduced  a  potential  function  and  associated  coordinate  system 
that  would  lead  to  the  separability  of  the  Hamilton-Jacobi  equations.  In  a  subse¬ 
quent  work  [Ref.  11],  Vinti  showed  that  the  equations  result  in  a  solution  reducible 
to  quadratures.  The  Vinti  potential  is  an  exact  expression  of  the  Earth’s  potential 
through  the  second  harmonic.  The  theory  provides  perturbed  coordinates,  not 
perturbed  elements.  In  his  second  order  analysis  [Ref.  12],  Kozai  criticized  this 
characteristic,  and  he  chose  not  to  use  the  Vinti  potential  since  it  would  require 
changing  the  definitions  of  the  conventional  orbital  elements. 
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Morrison  [Ref.  13]  showed  that  the  von  Zeipel  method  is  a  particular  case  of 
the  method  of  averaging,  and  Liu  [Ref.  14]  used  the  latter  technique  to  study  the 
combined  effects  of  air  drag  and  planet  oblateness  on  a  satellite  orbit.  The  method 
of  averaging,  unlike  von  Ziepel’s  method,  does  not  require  that  transformations 
be  canoi  ‘~al.  The  method  of  averaging  has  been  used  extensively  in  recent  years 
by  Lorell  and  Liu  (Ref.  15],  McClain  [Ref.  16],  and  Hoots  [Ref.  17].  However,  the 
validity  of  the  method  has  been  challenged;  most  notably  Taff  [Ref.  18]  doubts 
its  rigorous  foundation.  Arnold  [Ref.  19]  notes  that  the  principle  of  averaging  is  a 
vaguely  formulated  and  rigorously  untrue  assertion,  but  he  adds  that  sometimes 
such  assertions  are  fruitful  mathematical  sources.  It  should  be  also  noted  that 
the  solutions  obtained  in  Liu’s  analysis  [Ref.  14]  are  not  valid  for  near  circular 
equatorial  orbits  nor  at  the  critical  inclination,  while  those  obtained  by  Hoots  are 
valid  only  for  small  eccentricities  (0  <  e  <  .1),  and  they  are  invalid  at  inclinations 
near  0°  or  the  critical. 

Hamiltonian  mechanics  has  provided  a  rich  source  of  literature  in  orbit  the¬ 
ory;  however  its  practical  applicability  has  been  questioned  in  some  textbooks. 
Roy  [Ref.  20]  briefly  outlines  the  use  of  the  canonic  equations,  while  Taff  [Ref. 
18p.  322]  states  tnat  he  does  not  see  any  additional  practical  applications  provided. 
Baker  [Ref.  21]  chooses  not  to  represent  the  subject.  A  general  criticism  is  that 
the  process  of  generating  suitable  transformations  in  the  perturbation  procedure 
tends  to  make  the  coordinates  and  the  momenta  less  distinguishable  on  physical 
grounds  and  more  difficult  to  relate  to  the  set  of  natural  coordinates  which  were 
used  to  write  down  the  initial  set  of  differential  equations.  While  the  ultimate 
form  of  the  governing  equations  may  be  simple  to  solve,  there  remains  the  tedious 
task  of  obtaining  explicit  results  in  terms  of  physically  meaningful  coordinates  or 
elements. 
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C.  KOZAI’S  METHOD 

Prior  to  his  work  cited  above,  Kozai  [Ref.  22]  developed  a  method  for  find¬ 
ing  the  perturbations  on  the  orbital  elements  of  a  satellite  considering  only  the 
oblateness  of  the  Earth.  Kozai  developed  a  disturbing  function  based  on  the 
Earth’s  departure  from  a  sphere,  and  he  used  a  version  of  Lagrange’s  planetary 
equations  to  formulate  the  solution.  Kozai  used  the  standard  form  for  the  Earth’s 
potential  and  included  the  harmonics  J2,  Js,  and  J4.  Despite  the  use  of  the  higher 
harmonics,  the  theory  is  first  order.  Kozai  expressed  short-period  terms  in  J2,  the 
secular  in  J2,  J4,  and  J|,  and  the  long  period  terms  in  J2,  and  j4.  The  analytic 
expressions  are  developed  using  the  standard  orbital  elements. 

Kozai’s  work  is  cited  here  because,  due  to  its  simplicity,  the  method  has  be¬ 
come  very  popular  in  many  textbooks  and  handbooks  on  orbital  mechanics.  [Refs. 
20,  23,  24].  However,  TafF  cites  Kozai’s  method  as  an  example  of  misapplication 
of  perturbation  theory  [Ref.  18p.  332] .  Taff  challenges  the  assumptions  made  by 
Kozai  in  his  analysis,  and  he  points  out  that  the  method  is  invalid  at  the  critical 
inclination. 

As  was  stated  in  Chapter  I,  a  motivation  for  this  current  analysis  was  the 
purpose  of  finding  a  perturbation  method  that  would  lead  to  universal  solution.  All 
of  the  methods  discussed  thus  far  have  particular  problems  at  certain  inclinations 
or  eccentricities.  Some  of  the  problem  cases  have  been  resolved  by  unique  efforts 
(These  cases  are  discussed  in  Appendix  C.);  however,  one  should  question  the 
underlying  validity  of  any  perturbation  method  that  produces  singularities.  There 
has  been  no  satisfactory  way  found  for  avoiding  the  critical  inclination  singularity 
in  Kozai’s  method. 
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D.  THE  DIRECT  METHODS 


R.  E.  Roberson  [Ref.  25]  devised  an  approach  for  finding  the  qualitative 
and  approximate  quantitative  results  concerning  the  behavior  of  a  set  of  orbital 
elements  in  the  gravitational  field  of  an  oblate  planet.  The  motion  of  a  near 
satellite  around  the  planet  is  simply  described  to  first  order  by  introducing  a  frame 
of  reference  which  contains  a  mean  orbital  plane  having  a  constant  inclination. 
Both  the  reference  frame  and  the  orbit  plane  rotate  at  a  constant  angular  velocity. 
With  respect  to  this  doubly  moving  reference  frame,  the  motion  of  the  satellite 
differs  from  pure  elliptic  motions  by  only  periodic  perturbations. 

King-Hele  [Ref.  26]  advanced  the  approach  taken  by  Roberson.  He  intro¬ 
duced  a  non-rotating  reference  frame  with  a  orbit  plane  that  continually  rotates 
at  a  non-constant  rate  about  the  Earth’s  axis.  A  relation  is  found  between  the 
rotating  orbit  plane  and  the  angular  rate  of  travel  of  the  satellite.  The  equations 
of  motion  are  written  in  position  coordinates,  but  are  subsequently  rearranged 
in  terms  of  a  modified  set  of  orbital  elements.  The  inclination  of  the  rotating 
reference  plane  is  held  strictly  constant  in  the  analysis.  King-Hele  formulated 
the  problem  by  employing  a  power  series  expansion  in  terms  of  the  eccentricity; 
therefore,  the  method  is  limited  to  small  eccentricities.  The  final  solution  con¬ 
tains  an  incomplete  set  of  workable  elements  and  hence  the  method  gives  mostly 
a  qualitative  description  of  the  satellite’s  behavior  due  to  oblateness. 

King-Hele’s  analysis  was  the  inspiration  for  the  work  of  Brenner  and  Latta 
[Ref.  27j.  They  improved  King-Hele’s  original  analysis  by  abandoning  the  condi¬ 
tion  of  constant  orbit  inclination  and  by  retaining  the  eccentricity  in  closed  form 
expressions.  Brenner  and  Latta  limited  their  analysis  to  small  eccentricity  al¬ 
though  the  method  is  not  so  restricted.  They  obtained  an  approximate  first  order 
solution  and  demonstrated  that  the  method  is  valid  for  higher  order  analysis. 
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An  advantage  of  the  direct  method  is  that  one  may  use  ordinary  perturbation 
analysis  (the  Method  of  Strained  Coordinates)  to  solve  the  equations  of  motion.  In 
addition  the  chosen  set  of  orbital  elements  used  throughout  the  analysis  correspond 
closely  with  the  classical  elements;  therefore,  physical  interpretation  is  facilitated. 

The  disadvantage  in  using  the  method  is  shared  by  most  all  other  procedures 
that  must  use  a  perturbation  scheme.  The  process  requires  the  manipulation  of 
massive  algebraic  expressions.  However,  this  drawback  has  been  greatly  reduced 
by  the  introduction  recently  of  large  symbolic  mathematics  programs,  such  as 
MACSYMA,  that  handle  the  bookkeeping. 
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III.  THE  TWO-BODY  PROBLEM 

A.  INTRODUCTION 

The  central  problem  of  celestial  mechanics  is  the  two-body  (or  Kepler)  prob¬ 
lem.  Simply  stated,  the  problem  is  to  solve  for  the  motion  of  two  particles  inter¬ 
acting  through  their  mutual  gravitation  in  an  isolated  space.  A  solution  of  the  two 
body  problem  often  represents  physical  reality  in  an  acceptable  way.  For  instance, 
the  orbit  of  the  Earth  around  the  Sun  may  be  treated,  to  a  6rst  approximation,  as 
a  two-body  problem  because  the  influeni  e  of  perturbing  bodies,  the  Moon,  Jupiter, 
etc.,  are  small  compared  with  the  Sun's  gravitational  attraction.  Likewise,  more 
complex  problems  such  as  a  spacecraft  mission  to  Mars,  a  four  body  problem  - 
Earth,  Sun,  Mars,  spacecraft,  may  be  treated  by  breaking  the  flight  into  three 
two-body  problems.  Such  a  technique  is  used  in  the  patched  conic  method  where 
two-body  solutions  are  literally  patched  together. 

There  are  other  more  important  reasons  for  studying  the  two-body  problem. 
It  is  the  only  gravitational  problem  in  dynamics,  other  than  very  specialized  cases 
in  the  three-bodv  problem,  for  which  there  is  a  complete  and  general  solution, 
and  it  is  possible  to  gain  considerable  insight  into  more  general  phenomena  of 
motion  by  a  thorough  study  of  the  two-body  problem.  In  fact,  the  most  complete 
theories  of  celestial  motion  use  functions  appearing  in  the  solution  of  the  two- 
body  (elliptic  case)  problem  as  elementary  functions.  The  solution  is  central  to 
this  present  analysis  because  it  will  serve  as  a  starting  point  for  the  generating  of 
analytical  solutions  that  are  valid  to  higher  orders  of  accuracy.  These  solutions, 
called  general  perturbation  theories,  are  the  subject  of  Chaper  VI. 
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Li  this  chapter,  the  equations  for  the  two-body  problem  will  be  derived  and 
solved.  The  first  step  is  to  chose  a  coordinate  system  in  which  the  laws  of  Newton 
hold  (an  inertial  coordinate  system).  In  practice,  the  reference  frame  of  the  “fixed” 
stars  provides  a  very  good  approximation  to  an  inertial  reference  frame.  Next, 
following  the  method  of  Nelson  and  Loft  [Ref.28  pp. 82-84]  it  will  be  shown  that 
the  center  of  mass  of  two  bodies  in  this  coordinate  system  is  also  an  inertial 
reference  point.  Then  the  equations  of  motion  will  be  derived  and  solved  using 
the  polar  angle  6  as  the  independent  variable. 

B.  THE  DIFFERENTIAL  EQUATION 

Figure  3.1  shows  two  mass  centers  at  position  Tj  and  r2.  The  o/igin  O  is 
defined  to  be  an  inertial  reference  point.  The  distance  between  the  two  mass 
centers  is  r  where 

r  =  r2  -  n. 

Combining  Newton’s  third  law  of  motion  and  his  law  of  universal  gravitation  gives 

K 


Figure  3.1:  Two  Body  Problem 
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the  force  equations  for  the  two  bodies: 

d?Ti  Gmi  m2  r 


mj 


mj- 


dt2 
(Pti 


r3 

-Gmimj  r 


dt2  rs 

where  G  is  the  universal  gravitational  constant. 
Adding  (3.1)  and  (3.2)  and  integrating  results  in 


(3.1) 

(3.2) 


dri  dr2  .  . 

mi  — — (-  mj—  =  constant.  (3-3) 

dt  dt 

Now  the  center  of  mass  of  the  system  (barycenter)  in  Fig.  3.1  is  defined  as: 


(mi  +  m2)r0  =  miri  +  m2r2.  (3.4) 

By  combining  equations  (3.3)  and  (3.4)  the  following  result  is  obtained, 

dr0  mi  dri  m2  dr2 

— —  = - - — I — - —  =  constant. 

dt  m i  +  m2  dt  mi  +  m2  dt 

So  the  barycenter  moves  with  constant  velocity,  and  it  too  is  an  inertial  reference 
point.  Subtracting  (3.1)  from  (3.2)  results  in 


d2  r2  d2Ti  d?r  — G(mi  +  m2)r 


(3.5) 


dt 2  dt 2  dt 2  r3 

the  solution  of  which  gives  the  position  of  either  body  relative  to  the  other.  Now 
choose  the  barycenter  as  the  origin  and  define  the  position  of  m2  and  m2,  respec¬ 
tively,  as 

— m2 


r0i  =  rj  -  rc,  = 


mi  +  m2 


ro2  —  r2  —  r&  — 


mj 


mi  +  m2 


-r. 
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Substitution  of  the  above  equations  in  (3.1)  and  (3.2)  yield  separate  equations  for 
the  motion  of  each  body  relative  to  the  barycenter: 


<Pt0 

,  ^roi 

(Ptoi 

dt2 

j - it 

dt2 

dt2 

d?r0 

(PT02 

d2r02 

dt 2 

dt2 

dt2 

Equation  (3.5)  and  equations  (3.6)  are 
stant,  so  that 


—G  m\ 

(mx  -I-  mj)^!  f01 

(3.6) 

-Gm\ 

(mi  +  m2)2  r^2  r°2' 

identical  form,  differing  only  by  a  con- 


d?T  GMt 
+  — r-  =  0 


(3.7) 


dt2 

is  the  vector  differential  equation  of  motion  for  either  of  the  two  bodies,  r  is  the 
distance  to  the  other  body,  or  to  the  barycenter,  according  to  the  appropriate 
choice  of  GM.  For  the  problem  of  a  satellite  orbiting  the  earth,  the  mass  of  the 
satellite  can  be  neglected  in  comparison  with  the  mass  of  the  earth,  therefore  GM 
is  the  product  of  the  universal  constant  and  the  mass  of  the  earth. 


C.  THE  INTEGRATION  OF  THE  TWO-BODY  PROBLEM 

There  are  no  cross  products  involved  in  equation  (3.7);  therefore,  all  motion 
must  lie  in  the  plane  that  contains  r  and  The  scalar  components  of  acceleration 
are 


d*0  ^  2dr  d0 

dt 2  dt  dt 


(3.8) 


d2r  fd6\2  _-GM 
dt2  \  dt )  r2 


(3.9) 
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Writing  (3.8)  as  “  (f2af)  =  0  and  integrating  yields 


f  dt 


=  h  = 


constant 


(3.10) 


where  h  is  the  specific  angular  momentum. 

Equation  (3.10)  is  an  exact  integral  of  (3.8).  It  corresponds  to  Kepler’s 
empirical  law  of  constant  areal  velocity  which  states  that  the  area  swept  out  by 
the  radius  vector  of  a  planet  is  uniform  in  time. 

From  (3.10),  the  independent  variable  may  be  changed  from  t  to  6 ,  e.g., 


_  h_  d_ 
dt  r2  d6 

so  (3.9)  becomes 

1  d2r  2  fdr\2  1  GM  d 2  „  ,  1  GM 

r2d62  +  r*{d8j  +  r  ~  h2  ~d8^l'r’  +  r~  h2  ' 

Since  this  equation  is  linear  in  the  reciprocal  of  r,  it  may  be  written  as 


d?u  GM 

dd2^U~  h2 

which  has  the  solution 


GM 

u  =  — — — i-  v4cos(0  -  w)  (3-11) 

hr 

or  reintroducing  r,  (3.11)  becomes 

_ h2  /GM _ 

1  +  (Ah2 /GM)  cos (6  -  t v) 


which  may  be  written  as  the  polar  equation  of  a  conic  section 


r  = 


(3.12) 


P 

1  +  e  cos(0  -  u) 

so  that 

p  =  h2/GM  and  e  =  Ah*/GM. 

D.  ELLIPTIC  MOTION 

Equation  (3.12)  is  the  equation  of  a  conic  with  prime  focus  at  O.  The  conic 
has  a  semi-latus  rectum  h2/GM,  an  eccentricity  e,  and  a  semi-major  axis  a  that 
makes  an  angle  f  —  6  —  u  with  the  horizontal  axis  (Figure  3.2).  The  extreme 
endpoints  of  the  major  axis  of  the  orbit  are  referred  to  as  apsides  or  apses.  The 
point  nearest  the  prime  focus  is  called  perigee  and  is  given  by  6  =  u>.  The  point 
farthest  from  the  prime  focus  is  given  by  0  =  w  +  180°  and  is  called  apogee.  The 
angle  u,  “argument  of  perigee”,  will  be  discussed  later  in  this  chapter. 

The  energy  of  the  satellite  in  the  orbit  is  conserved  and  is  equal  to 

E  =  m,(i)2/ 2  -  GM/r)  =  m,C 


Figure  3.2:  The  Elliptic  Orbit 
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where  C  is  the  total  energy  of  the  satellite,  t>*/2  the  kinetic  energy  and  —GM/r 
the  potential  energy  of  the  satellite,  all  per  unit  mass  (m#).  In  general  C  is  equal 
to  —GM/2a,  so  that  the  satellite’s  energy  depends  only  on  the  semi-major  axis. 
From  this  relationship  it  is  easily  shown  that  the  satellite’s  velocity  is  given  by 


v*  =  GM(2/r  —  1/a). 

Then  the  relationships  rp  =  a(l  —  e)  and  ra  =  a(l  -f-  e)  result  in 

j  GM{  1  +  e) 
v’~  o(l  -  e) 

and 

2  GM{  1  -  e) 

Va  ~  a{ 1  +  e) 

so  that  the  velocity  is  a  maximum  at  perigee  and  a  minimum  at  apogee. 

The  area  of  the  ellipse  is  na2^/(l  —  e2),  and  the  rate  of  description  of  area  is 
~~  =  Since  h 2  =  GMa(l  -  e2),  the  orbital  period  may  be  written  as 

™  2tt  s/, 

T  =  -y=a3/2. 

\[GM 

By  defining  the  mean  motion  n  as  T  =  so  that  n2a3  =  GMt  one  may  proceed 
to  derive  an  expression  for  position  versus  time  in  the  elliptic  orbit. 

The  orbital  ellipse  APB  with  center  at  C  touches  at  perigee,  A,  and  at 
apogee,  B ,  a  concentric  circle  also  centered  at  C  which  has  as  a  radius  the  semi¬ 
major  axis  of  the  ellipse.  The  circle  C  is  known  as  the  auxiliary  circle,  and  is 
geometrically  related  to  the  ellipse  by  the  relation 
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where  t  is  the  eccentricity,  P  any  point  on  the  ellipse,  N  the  foot  of  the  perpendic¬ 
ular  through  P  upon  AB ,  and  P'  the  intersection  of  this  perpendicular  with  the 
circle  C.  The  angle  ACP1  =  E  is  known  as  the  eccentric  anomaly  (Figure  3.3). 

Let  r  be  the  time  of  perigee  passage  and  t  the  time,  then  t  —  t  is  the  time 
since  perigee  passage.  The  quantity, 

M  =  n(t  —  t ),  (3.13) 

is  called  the  mean  anomaly.  Using  the  geometry  of  Figure  3.3,  one  may  derive  the 
following  relationship  between  the  anomalies: 

M  =  E  -  esin  E,  (3-14) 

which  is  known  as  Kepler’s  equation.  If  the  position  of  a  satellite  in  a  fixed  orbit 
relative  to  the  earth  is  desired  at  some  specified  time  t,  then  equation  (3.13)  gives 
M  and  (3.14)  gives  E  .  The  distance  to  the  satellite  is  found  by  the  relationship 

r  =  g(1  -  t  cos  E ). 

The  angular  position  of  the  satellite  is  defined  by  the  true  anomaly,  8  —  vj  =  f ,  or 
the  angular  distance  from  perigee: 


Figure  3.3:  The  Eccentric  Anomaly 
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tan  (//2)  =  ^  — -  J  tan-. 

In  actual  practice,  the  second  step  in  the  above  process,  solving  (3.14)  is  a  bit 
more  involved  since  a  closed  form  solution  to  Kepler’s  equation  does  not  exist. 
However,  dozens  of  methods  of  successive  approximations  have  been  devised.  For 
instance,  by  use  of  a  numerical  method  M  can  be  calculated  for  a  few  values  of  E 
and  then  the  correct  E  can  be  found  by  inverse  interpolation. 

E.  CONSTANTS  OF  THE  MOTION 

In  section  B,  the  original  equations  of  motion  (3.1)  and  (3.2)  were  reduced 
to  (3.7).  Thus  the  problem  was  reduced  from  one  of  three  second  order  differential 
equations  requiring  twelve  constants  to  one  of  three  second  order  equations  with  six 
constants.  A  discussion  of  the  constants  of  motion,  some  of  which  were  introduced 
during  the  solution  of  the  two-body  equation,  is  the  purpose  of  this  section.  The 
six  constants  may  be  written  in  a  variety  of  forms,  a  choice  among  forms  is  usually 
made  with  the  purpose  of  simplifying  the  problem. 

Equation  (3.7)  was  solved  by  the  classic  technique  of  changing  the  indepen¬ 
dent  variable  from  t  to  6.  But  (3.7)  in  its  present  form  is  integrable;  that  is, 
there  exists  sufficient  time  independent  first  integrals,  or  functions  that  are  con¬ 
stant  along  the  motion,  to  specify  each  orbit.  The  first  of  these  is  the  angular 
momentum.  Cross-multiplying  (3.7)  by  r  results  in 


or 


d2r  GM 
r  x  —  +  r  x  — ,-r  =  0 


dt2 


d2r 

t  x  —  =  0. 
dt 2 
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And  since 


then 


By  integration 


d  (  dr\  dr  dr  <Pr 


|(r  X  v)  =  0. 


rxv  =  h 


where  h,  the  angular  momentum,  provides  three  constants  of  motion.  Similarly, 
by  cross-multiplying  (3.7)  by  h,  one  may  obtain  the  Lenz  vector 


dr  h  r  /  de  \ 

®  ~  dtX  GM~  r  (di’=0j 

where  e  is  a  vector  along  the  major  axis  of  the  orbit  pointing  toward  the  position  of 
closest  approach  or  perigee  (|e|  =  e).  The  vector  e  provides  only  two  independent 
constants  since  h  and  e  are  perpendicular  vectors  (h  •  e  =  0),  and  one  remaining 
constant  is  required. 

The  vector  integrals  h  and  e  specify  that  the  orbit  will  lie  in  the  plane  per¬ 
pendicular  to  h  and  have  a  shape  determined  by  e.  The  classical  orbital  elements: 
a  (semi-major  axis),  e  (eccentricity),  i  (inclination),  fi  (longitude  of  the  ascending 
node),  and  u  (argument  of  perigee),  may  be  derived  from  these  vectors. 

From  (3.12)  the  equation  for  r  may  be  written  as 


h*/GM 

1  +  e  cos  / 


where  /  is  the  angle  between  e  and  r  (the  true  anomaly). 
Restricting  e  to  the  elliptic  case  results  in 
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fl  GM{  1  +  e)  ~  T  ~  GM{  1  -  e)  f* 

so  that  the  major  axis  of  the  ellipse  is 

o  2hi 

2o_r,  +  r’  =  GM( ITpj- 

Knowing  t  and  a  gives  the  shape  and  size  of  the  orbit.  The  orbit  now  can  be 
oriented  in  a  coordinate  system.  Reference  is  made  to  Figure  3.4  and  the  two 
angles  t  and  fl.  t  is  the  inclination  of  the  orbit  plane  defined  as  the  angle  between 
the  equatorial  plane  and  the  orbit  plane.  Since  this  is  the  same  angle  as  the  angle 
between  the  z  axis  (k  unit  vector)  and  the  angular  momentum  vector  h,  t  may  be 
found  by 


Figure  3.4:  Constants  of  the  Orbit 
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The  angle  0,  the  longitude  of  the  ascending  node,  is  the  angle  in  the  equato¬ 
rial  plane,  between  the  i  unit  vector  and  the  point  N  where  the  satellite  crosses 
through  the  equatorial  plane  in  a  northerly  direction  measured  counterclockwise 
when  viewed  from  the  north  side  of  the  equatorial  plane.  The  vector  n  lies  along 
NO  such  that 

n  =  k  x  h. 

Therefore  fl  may  be  found  by 

cos  n  =  — . 
n 

As  was  stated  above,  the  vector  e,  points  toward  perigee.  The  angle  u>  (the 
argument  of  perigee)  measures  the  distance  between  NO  and  perigee  and  can  be 
found  by 

n  ■  e 

COS  U)  =  - . 

ne 

With  the  constants  a,  e,t,  fi,w  specified,  the  orbit  is  defined  in  the  coordinate 
system.  The  remaining  task,  to  locate  the  satellite  in  the  orbit  at  any  time, 
requires  one  more  constant. 

The  final  constant  of  motion  is  given  by  the  relationship  between  the  mag¬ 
nitude  of  the  angular  momentum  and  the  true  anomaly  (/): 

[  hdt  =  /  r7df. 

JT  •'0 

By  a  change  of  variables  from  /  to  the  eccentric  anomaly  E,  the  above  equation 
may  be  easily  integrated.  The  result  by  this  analytical  method  is  identical  to 
equation  (3.14),  Kepler’s  equation,  which  was  previously  derived  by  geometry. 
The  constant  r,  time  of  perigee  passage,  is  the  final  constant  of  integration. 
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F.  SUMMARY 


This  completes  the  analysis  of  the  two-body  problem.  It  has  been  shown  that 
a  combination  of  six  constants  will  strictly  define  the  motion  of  a  satellite  under 
the  influence  of  a  central  gravitational  force.  The  six  that  were  chosen,  referred 
to  now  as  the  orbital  elements,  can  be  used  to  find  other  constants  of  the  motion 
including  r  and  v,  the  distance  and  velocity  vectors  of  the  satellite. 
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IV.  SATELLITE  PERTURBATIONS 


A.  INTRODUCTION 

As  was  demonstrated  in  the  last  chapter,  the  classical  two-body  problem  has 
solutions  that  can  be  written  in  closed  form  when  the  polar  angle  (or  the  eccen¬ 
tric  or  true  anomaly)  is  used  as  the  independent  variable.  If  an  additional  force 
acting  on  either  of  the  two  bodies  is  introduced,  the  resulting  equations  of  motion 
usually  no  longer  have  closed-form  solutions.  When  the  magnitude  of  such  a  force 
is  small  compared  to  the  central  gravity  term,  the  force  is  called  a  perturbation. 
The  resulting  orbit  experiences  a  small  departure  from  the  Keplerian  orbit,  at 
least  initially.  These  departures  are  also  called  perturbations.  Under  certain  cir¬ 
cumstances,  it  is  possible  to  make  analytic  approximations  to  the  effects  of  the 
perturbing  forces,  though  a  precise  solution  cannot  be  obtained.  Generally,  the 
methods  consist  in  determining  the  exact  equations  of  motion  and  then  assuming 
that  their  solutions  do  not  depart  appreciably  from  the  case  of  no  disturbing  force. 
Then  only  an  indication  of  the  actual  motion  of  the  body  can  be  obtained.  Pre¬ 
cise  solutions  can  be  found  for  specific  initial  conditions  by  numerical  integration 
techniques,  but  these  solutions  give  little  insight  into  the  dependence  of  the  mo¬ 
tion  on  the  parameters  of  the  disturbing  force.  In  some  cases,  the  approximations 
obtained  with  analytic  methods  may  exceed  the  precision  of  numerical  methods, 
especially  if  the  prediction  is  required  for  a  long  period  of  time  and  there  is  a  clear 
dominance  of  one  particular  force. 


In  the  case  of  a  close  satellite  about  a  non-spherical  planet,  a  potential  func 
tion  V  can  be  formed  such  that 


V  =  Vo +  R 

where  Vo  is  the  potential  function  due  to  the  two-body  problem  and  R  the  dis¬ 
turbing  function  that  is  at  least  an  order  of  magnitude  less  than  Vo-  Many  general 
perturbation  theories  make  use  of  the  fact  that  the  two-body  orbit  of  the  body  due 
to  Vq  changes  slowly  due  to  /?,  and  they  attempt  to  obtain  analytical  expressions 
for  the  changes  in  the  orbital  elements  due  to  R  within  a  specified  time  interval. 
If  the  elements  of  an  elliptical  orbit  are  Oo,  e0,  t0,  n0,  tv0,  and  t0,  the  ellipse  with 
these  elements  is  referred  to  as  the  osculating  elements  at  to-  The  velocity  of  the 
disturbed  satellite  at  this  time  in  its  osculating  ellipse  is  equal  to  its  velocity  in 
the  actual  orbit. 

At  a  future  time  tx,  the  elements  will  change  due  to  the  presence  of  R. 
For  instance,  Oq  will  be  changed  to  aj,  and  the  quantities  (gj  —  o0),(ei  -  e0). 
etc.  are  called  perturbations  in  the  elements.  Corresponding  to  the  perturbations 
in  the  elements  are  perturbations  in  the  coordinates  and  velocity  components. 
There  are,  however,  at  least  two  reasons  for  using  orbital  elements  rather  than 
coordinates  to  describe  the  motion  of  the  satellite.  First,  the  elements  do  not 
exhibit  a  variability  of  anomalistic  motion  that  the  coordinates  do,  therefore  any 
variation  can  be  attributed  directly  to  the  perturbing  forces.  Second,  the  elements 
give  a  clearer  geometric  picture  than  do  the  coordinates;  hence  the  effect  of  the 
perturbation  on  the  orbit  can  be  seen  immediately. 

There  are  various  kinds  of  disturbances  that  an  orbit  can  experience,  the 
severity  of  each  is  usually  due  to  the  altitude  of  the  satellite.  It  is  the  purpose  of 
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this  chapter  to  give  a  qualitative  description  of  the  most  important  disturbances, 
and  to  relate  the  relative  magnitudes  of  each. 

B.  THE  EARTH’S  GRAVITATIONAL  FIELD 

The  two-body  problem  assumes  that  the  earth  is  a  sphere,  while  in  reality  the 
earth  is  flattened  somewhat  at  the  poles  and  bulges  correspondingly  at  the  equator. 
Such  a  shape  is  called  an  oblate  spheroid.  In  the  science  of  geodesy,  it  has  been 
useful  to  define  a  reference  ellipsoid  as  a  mathematical  surface  which  is  an  idealized 
approximation  to  the  earth’s  actual  surface.  The  study  of  satellite  orbits  has 
established  a  flattening  of  the  terrestrial  ellipsoid  as  1/298.24,  which  corresponds 
to  a  difference  between  the  equatorial  and  polar  radii  of  21.4  kilometers. 

Another  surface  commonly  used  in  geodesy  is  the  geoid,  which  is  the  equipo- 
tential  surface  that  coincides  on  the  average  with  mean  sea  level  in  the  oceans.  The 
geoid  is  everywhere  perpendicular  to  a  plumb-line  since  gravity  is  always  normal  to 
its  surface.  Before  the  advent  of  artificial  satellites,  it  was  generally  accepted  that 
except  for  relatively  small  variations  resulting  from  the  presence  of  mountains  or 
deep  valleys,  the  geoid  could  be  regarded  as  approximately  an  ellipsoid.  Now,  the 
surface  of  constant  gravitation  can  be  more  accurately  portrayed  by  representing 
the  potential  as  a  series  of  quantities  known  as  spherical  harmonics,  each  of  which 
makes  a  contribution,  positive,  negative,  or  zero,  to  the  total.  The  contribution 
of  any  harmonic  is  determined  by  a  factor,  usually  represented  by  the  symbol  J 
and  commonly  referred  to  as  the  value  of  that  harmonic.  These  J  values  for  a 
planet’s  gravitational  field  can  be  determined  from  observations  of  a  satellite  orbit 
and  they  can  be  related  to  the  shape  of  the  geoid. 

A  large  number  of  harmonics  may  be  required  to  precisely  represent  a  planet’s 
gravitational  field,  but  in  practice  the  higher  harmonics  make  such  a  small  contribution 
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that  they  can  be  neglected,  at  least  to  a  first  approximation.  The  zonal  harmonic 
J0  expresses  the  overall  size  of  the  geoid,  while  Jls  the  first  degree  harmonic  de¬ 
termines  the  center  point  of  the  geoid  in  the  north-south  direction.  The  other 
harmonics  represent  deviations  from  the  spherical  shape  as  shown  by  Figure  4.1. 
It  is  seen  that  the  contributions  from  the  even  harmonics  are  symmetrical  about 
the  equator,  while  the  odd  harmonics  corresponds  to  contributions  that  are  asym¬ 
metrical.  The  degree  of  the  harmonic  gives  the  number  of  undulations  in  the  shape 
of  the  surface. 


Figure  4.1:  Qualitative  Representation  Of  The  Harmonics  Of  The  Geoid 

The  results  given  thus  far  consider  only  the  zonal  harmonics,  which  are  in¬ 
dependent  of  longitude.  The  tesseral  harmonics  give  east-west  deviations  from 
symmetry.  Satellite  observations  of  the  tesseral  harmonics  have  led  to  the  con¬ 
clusion  that  the  equator  of  the  earth’s  geoid  is  slightly  elliptical,  the  difference 
between  the  longest  and  shortest  axes  being  about  400  meters. 

The  following  tables  from  Kozai  [Ref.  29j  gives  representative  values  for  the 
coefficients  of  the  earth’s  zonal  and  tesseral  harmonics. 
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TABLE  4.1:  ZONAL  HARMONICS 


n 

J„  x  106 

n 

Jn 

x  10' 

— 

2 

1082.48 

± 

.04 

6 

.39 

± 

.09 

3 

-2.57 

± 

.01 

7 

-.47 

± 

.02 

4 

-1.84 

± 

.09 

8 

-.02 

± 

.07 

5 

-.06 

± 

.02 

9 

.11 

± 

.03 

TABLE  4.2:  TESSERAL  HARMONICS 


n 

m 

•c 

x  106 

2 

2 

2.32 

± 

.30 

—37°. 5 

± 

5°.  6 

3 

1 

3.95 

± 

.36 

22°.0 

± 

11°.0 

3 

2 

.41 

± 

.36 

31°. 0 

± 

14°. 0 

3 

3 

1.91 

± 

.29 

51°. 3 

± 

2°.9 

4 

1 

2.64 

± 

.44 

163°. 5 

± 

6°  .5 

4 

2 

.17 

± 

.06 

54°.0 

± 

11°.0 

4 

3 

.05 

± 

.04 

— 13°.0 

± 

19°  .0 

C.  ATMOSPHERIC  DRAG 

Atmospheric  drag  is  the  most  complex  and  the  most  difficult  of  the  artificial 
satellite  perturbations.  The  complexity  arises  from  the  fact  that  an  exact  force 
law  is  not  known  and  the  atmosphere  is  variable.  This  variability  results  from 
the  fact  that  the  atmosphere  is  not  spherically  symmetric  and  that  lunisolar  tides, 
diurnal  heating,  strength  of  the  solar  wind,  etc.,  all  effect  atmospheric  density. 
In  addition,  the  atmosphere  is  moving  with  the  rotating  earth.  Therefore  this 
perturbation  is  difficult  to  model  analytically  or  numerically. 

For  most  Earth  satellites,  atmospheric  drag  removes  the  satellite’s  energy, 
thereby  causing  it  to  drop  in  altitude  and  thus  increase  its  speed.  The  increased 
speed  and  the  lower  altitude  increases  the  drag  force  with  the  result  that  the 
satellite  spirals  into  Earth. 
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Satellite  drag  is  measured  along  the  same  direction  as  the  velocity  vector  of 
the  satellite.  The  formula  used  is 

F d  =  ~^CDAp\v  -  v.|(v.  -  v) 

where  v  is  the  satellite’s  inertial  velocity,  va  is  the  inertial  velocity  of  the  atmo¬ 
sphere,  p  is  the  atmospheric  density,  A  is  the  effective  cross-sectional  area  of  the 
satellite,  and  Co  is  the  drag  coefficient,  Co  «  2. 

D.  SOLAR  RADIATION  PRESSURE 

In  contrast  to  drag  which  is  a  tangential  force,  solar  radiation  pressure  is 
radial.  By  solar  radiation  on  a  satellite  is  meant  the  net  acceleration  resulting 
from  the  interaction  (i.e.,  absorption,  reflection,  or  diffusion)  of  the  sunlight  with 
the  surface  of  the  satellite.  In  general,  the  illumination  of  a  low  altitude  satellite 
orbit  will  be  non-uniform  due  to  the  Earth’s  shadow,  albedo,  and  atmospheric 
absorption.  The  magnitude  of  the  solar  radiation  pressure  is  given  by 

AP 

AirmcD2 

where  A  is  the  effective  cross-sectional  area  of  the  satellite,  P  is  the  total  radiated 
solar  power,  m  is  the  satellite’s  mass,  c  is  the  speed  of  light,  and  D  is  the  satellite- 
Sun  distance. 

It  can  be  seen  that  the  perturbing  effect  of  solar  radiation  particularly  effects 
satellites  with  a  large  area  to  mass  ratio.  Such  was  the  case  of  the  satellite  Echo  I, 
which  when  inflated  was  a  sphere  with  a  total  exposed  area  of  about  31,000  square 
feet  and  weighed  only  135  pounds.  Its  initial  orbit  was  approximately  circular  with 
an  altitude  of  1000  miles,  but  the  pressure  of  solar  radiation  brought  down  the 
perigee  to  600  miles  at  times. 


A.  OTHER  PERTURBING  FORCES 


Gravitational  perturbations  due  to  other  celestial  bodies  (the  Moon,  the  Sun, 
and  the  planets)  are  caused  by  the  differences  between  the  force  on  the  Earth  and 
that  on  the  satellite  or  the  “tidal”  force.  Perturbing  forces  from  the  Sun  and  Moon 
become  significant  as  the  altitude  of  the  satellite  increases.  Planetary  perturba¬ 
tions  are  very  small,  with  Venus  and  Jupiter  providing  the  largest  contributions. 

Since  classical  mechanics  is  used  in  this  analysis,  relativistic  corrections  may 
be  regarded  as  a  small  perturbation  to  the  motion.  Other  minor  perturbing  forces 
include  atmospheric  lift  and  electromagnetic  forces. 

The  following  table  from  Milani  [Ref  30]  gives  the  magnitudes  of  various 
perturbations  on  three  different  satellites: 

TABLE  4.3:  ACCELERATIONS  ON  SPACECRAFT  AT  VARIOUS 

ALTITUDES  (cm/sec2) 


Cause 

Geosynchronous 
satellite 
a  =  42,160  KM 
A/M  =  .1  cm2/g 

LAGEOS 

12,270 

.007 

SEASAT 

7,100 

.2 

Earth’s 

monopole 

2.2  x  101 

2.8  x  102 

7.9  x  102 

Earth’s 

oblateness 

7.4  x  10~4 

1.0  x  10"1 

9.3  x  lO^1 

Higher 

harmonics 

4.5  x  10-10 

8.8  x  10~6 

7.0  x  10~4 

Moon 

7.3  x  10"4 

2.1  x  10“4 

1.3  x  10~4 

Sun 

3.3  x  10~4 

9.6  x  10  5 

5.6  x  10~5 

Venus 

4.3  x  10-8 

1.3  x  10"8 

7.3  x  10~9 

General 

Relativity 

2.3  x  10'9 

9.3  x  10"8 

4.9  x  10"7 

Air  Drag 

0  (?) 

3.0  x  10"10 

2.0  x  10~5 

So'ar 

Radiation 

4.6  x  10"6 

3.2  x  10"7 

9.2  x  10~6 

Earth’s 

albedo 

4.2  x  10-8 

3.4  x  10“8 

3.0  x  10~6 

V.  DEVELOPMENT  OF  THE  EQUATIONS  OF 

MOTION 


A.  PRELIMINARIES 

1.  Overview 

In  this  chapter,  the  groundwork  will  be  laid  for  solving  the  equations  of 
motion  for  a  satellite  about  an  oblate  planet.  To  begin,  a  discussion  of  the  assump¬ 
tions  made  in  the  analysis  will  be  given.  Second,  the  special  coordinate  system, 
which  was  introduced  by  King-Hele  [Ref.  26]  and  refined  by  Brenner  and  Latta 
[Ref.  27]  will  be  developed,  followed  by  a  derivation  of  the  relationships  among 
the  astronomical  angles  of  the  coordinate  system  and  the  angles  of  a  spherical 
coordinate  system.  Next,  an  expansion  for  the  potential  of  a  planet  modeled  as 
an  oblate  spheroid  will  be  derived.  Finally,  the  equations  will  be  transformed  so 
that  they  can  be  solved  by  an  ordinary  perturbation  method. 

2.  Assumptions  and  Limitations 

The  basis  for  the  assumptions  made  in  this  analysis  is  the  requirement 
that  a  solution  to  the  equations  of  motion  for  a  satellite  orbiting  a  planet  be 
accurate  to  within  a  relative  error  <  10-6.  Therefore,  all  perturbative  forces 
that  are  of  magnitude  10-6  or  smaller  compared  to  1  may  be  neglected.  This 
requirement  then  allows  one  to  model  the  earth  as  an  oblate  spheroid  with  an 
axially  symmetric  gravitational  potential.  This  is  to  say  that  all  zonal  harmonics 
except  th-  second  in  the  expanded  potential  may  be  neglected.  In  addition,  all 
coefficients  of  the  tesseral  harmonics  are  small  enough  to  neglect. 
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The  assumption  of  an  axially  symmetric  potential  is  a  reasonable  as¬ 
sumption  if  the  earth  is  used  as  an  example  of  an  oblate  planet.  As  was  suggested 
in  Chapter  IV,  if  one  normalizes  the  expanded  geopotential  of  the  earth  such  that 
the  contribution  from  the  central  force  (1/r*)  is  1,  then  the  contribution  from  the 
zonal  harmonic  for  the  oblateness  contribution  is  «  10-s,  and  all  other  harmonics 
(zonal  and  tesseral)  are  <  10-6. 

The  gravitational  effects  of  the  Sun  and  Moon  are  also  neglected.  This 
assumption  will  remain  valid  ov't  t.o  a  distance  of  a+  least  6, (XX)  kilometers  atv  .  <. 
the  Earth,  where  the  oblateness  perturbation  is  still  about  1000  times  larger  than 
the  attraction  of  the  Sun.  However,  at  altitudes  near  geosynchronous  (35,800  km) 
the  perturbative  forces  of  the  Sun  and  Moon  are  nearly  equal  to  that  of  oblateness, 
and  therefore  would  have  to  be  considered. 

The  most  important  limitation  to  the  analysis  is  the  neglect  of  atmo¬ 
spheric  drag.  For  high  altitude  satellites  with  small  eccentricity,  the  neglect  of  air 
drag  is  unimportant;  however,  for  very  low  altitude  satellites  the  effect  of  air  drag 
would  begin  to  dominate  the  oblateness  force  after  a  number  of  revolutions.  Also, 
high  altitude  satellites  with  highly  eccentric  orbits  would  be  greatly  affected  by 
drag  as  they  pass  through  perigee.  For  these  particular  cases,  the  desired  accuracy 
for  1000  revolutions  would  not  be  achievable:  however,  in  many  cases  the  solution 
could  be  accurate  for  a  few  orbits.  As  was  discussed  in  Chapter  IV,  air  drag  re¬ 
mains  the  most  difficult  perturbation  to  model.  Since  an  accurate  geopotential 
model  allows  for  a  better  determination  of  fluctuations  in  drag  of  the  atmosphere 
through  which  the  satellite  is  moving,  the  analysis  conducted  here  is  valuable  even 
when  air  drag  becomes  the  dominating  factor. 

Also  neglected  in  this  analysis,  are  the  remaining  perturbative  effects 
mentioned  in  Chapter  IV.  Solar  radiation  pressure,  the  Earth's  albedo  radiation 
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pressure,  and  the  general  relativistic  correction  to  Newton’s  law  of  gravity  are  all 
neglected  since  their  relative  contributions  are  small. 

The  final  limitation  to  the  analysis  is  the  assumption  that  the  differential 
equations  which  describe  the  motion  of  the  satellite  have  a  solution  asymptotic  to 
the  form  specified  in  the  analysis.  As  long  as  the  true  solution  does  not  depart 
from  the  form  specified,  the  analysis  is  valid. 

3.  Order  of  Magnitudes 

Throughout  this  analysis,  reference  will  he  made  to  the  relative  mag¬ 
nitudes  of  the  perturbing  forces;  hence  it  is  important  to  establish  a  convention 
for  orders  of  magnitudes.  The  approach  established  by  King-Hele  [Ref  26]  will  be 
followed  here  with  the  exception  that  a  small  eccentricity  is  not  assumed.  Here, 
J  «  J-i,  the  coefficient  of  the  second  zonal  harmonic  in  the  potential  is  used  as  the 
basic  small  unit.  Mathematically, 

f(JJ)  =  0{Jm6n) 

means  that  there  exists  for  all 

0  <  J  <  10'3  and  0  <  6  <  (2tt)103 
a  K  independent  of  J  and  6  such  that 

|/i  <  KJm6n 

The  following  terms  are  defined: 

Zero  order  =  O(l)  ~  1.0 

First  order  h  0(J),0(J20),0(J362), . . .  ~  10-3 

Second  order  =  0(J2),0(J36),0(J*62), . . .  ~  10“6 


B.  THE  COORDINATE  SYSTEM 


Figure  5.1  is  an  inertial  reference  system  of  spherical  coordinates  with  the 
prime  direction  pointing  toward  the  vernal  equinox  at  epoch  1900.0.  The  equato¬ 
rial  plane,  latitude  (6),  and  longitude  (<£),  as  shown  in  the  figure  have  their  usual 
meanings. 


Figure  5.1:  The  Reference  System 


As  was  demonstrated  in  Chapter  III,  the  path  of  a  satellite  of  a  strictly  spherical 
planet  lies  entirely  in  a  fixed  plane,  and  the  motion  of  the  satellite  is  described  by 
the  solution  to  the  two-body  problem.  With  the  angular  momentum  vector  fixed 
in  space,  the  intersection  of  the  orbit  plane  and  the  equatorial  plane  describe  the 
fixed  angle  fl  measured  from  the  prime  direction  to  the  intersection  (or  node). 
The  longitude  of  the  ascending  node  (H),  the  inclination  (»),  and  the  argument  of 
perigee  (u;),  fix  the  orbit  plane  in  the  coordinate  system. 
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The  effect  of  the  oblateness  perturbation  is,  in  general,  to  move  the  satellite 
out  of  the  original  two-body  orbital  plane.  One  specific  effect  is  to  rotate  the 
orbit  plane  about  the  planet  in  the  opposite  direction  to  the  satellite’s  motion. 
The  physics  of  this  motion  is  easily  demonstrated  in  Figure  5.2.  Oblateness  is 
represented  by  an  additional  mass  about  the  equator  of  the  planet.  The  radius 
vector  r  and  the  satellite  5  lie  in  a  plane  which  is  named  here  the  “reference 
plane” .  If  there  were  no  equatorial  bulge,  the  direction  of  the  gravitational  force 
would  coincide  with  the  radius  vector  end  the  angular  momentum  vector  would 
remain  in  a  constant  direction  normal  to  the  plane.  The  equatorial  bulge;  however, 
adds  a  component  of  force  that  does  not  lie  along  r.  This  additional  force  adds  a 
torque  f  =  rxF.  The  direction  of  r  is  into  the  paper  at  5  for  a  prograde  orbit. 
Therefore  as  the  satellite  orbits,  the  angular  momentum  vector  rotates  about  the 
Z  axis.  The  reference  plane  also  rotates,  and  the  rotation  is  measured  by  the 
change  in  ft. 


Poior 

axis 


Figure  5.2:  Rotation  of  the  Reference  Plane 


The  line  NON'  in  Figure  5.1  is  referred  to  as  the  “line  of  nodes”,  and  the 
ascending  node  describes  the  point  where  the  satellite  enters  the  northern  hemi- 
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sphere.  For  the  two-body  problem  with  a  non-rotating  plane,  the  line  of  nodes  is 
a  fixed  reference  line.  Now,  for  the  rotating  plane,  when  the  satellite  is  actually 
at  a  node,  the  line  of  nodes  passes  through  the  satellite,  but  at  other  times  the 
definition  “line  of  nodes”  is  an  arbitrary  one,  and  the  angle  O  is  simply  given  as 
a  function  of  time. 

Another  motion  of  the  reference  plane  that  is  a  result  of  the  oblateness 
perturbation  is  an  in-plane  rotation  called  precession.  As  was  demonstrated  in 
Chapter  III,  a  particle  executing  bounded,  non-circular  motion  in  a  central  force 
field  will  always  have  a  radial  distance  from  the  force  center  that  is  bounded  by 
rmax.  <  r  <  rm,„.  That  is,  r  is  bounded  by  the  apsidal  distance.  Such  an  orbit  is 
called  a  closed  orbit  and  it  is  characterized  by  the  fact  that  all  apsidal  angles  are 
equal.  For  instance,  the  apsidal  angle  is  rr  for  elliptic  motion.  But  if  the  radial 
dependence  deviates  slightly  from  1  /r*,  then  the  apsides  will  precess  or  rotate 
slowly  in  the  plane  of  motion  (Figure  5.3).  This  motion  is  analogous  to  the  slow 
rotation  of  elliptic  motion  of  a  two  dimensional  harmonic  oscillator  whose  natural 
frequencies  for  each  dimension  are  almost  equal.  The  rotation  of  the  line  of  apsides 
is  also  known  as  the  precession  of  perigee  since  the  value  for  6  for  which  r  is  a 
minimum  varies  in  time  as  the  apsides  rotate. 


Figure  5.3:  Rotation  of  the  Line  of  Apsides 
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For  a  satellite  orbiting  an  oblate  planet  the  rotation  of  the  line  of  apsides  is 
opposite  the  direction  of  satellite  motion  for  inclinations  less  than  sin-1  4/5  and 
in  the  same  direction  as  the  satellite  for  inclinations  greater  than  sin-1  ^4/5.  If 
the  inclination  is  exactly  sin-1  yji/ 5  then  there  is  no  rotation  of  the  line  of  apsides, 
and  perigee  remains  constant.  This  inclination  is  known  in  Astrodynamics  as  the 
“critical  inclination” .  A  detailed  discussion  of  the  critical  inclination  is  contained 
in  Appendix  C. 

A  final  motion  of  the  reference  plane  caused  by  the  oblateness  perturbation 
is  the  periodic  variation  of  the  inclination  about  the  initial  inclination  i0. 


C.  ANGLE  RELATIONS 

As  a  preliminary  to  deriving  the  equations  of  motion,  it  is  necessary  that 
the  relationships  among  the  spherical  coordinates  and  the  angles  describing  the 
reference  plane  be  established. 

Reference  is  made  to  Figure  5.4,  where  i,i',j,  and  j'  define  the  equatorial 
plane.  Let  a,b,  and  c  denote  three  unit  vectors,  where  b  and  c  are  in  the  reference 
plane,  c  points  to  the  initial  point  6  =  7t/2  in  the  reference  plane  where  6  is 
measured  from  the  line  of  nodes,  b  points  to  the  position  of  the  satellite,  S.  a  is 
perpendicular  to  the  reference  plane  and  therefore  points  in  the  same  direction  as 
the  angular  momentum  vector  h.  a  and  c  are  both  perpendicular  to  the  line  of 
nodes.  Then 


b  =  cos  <t>  cos  6  i  +  sin  <fr  cos  S  j  4-  sin  S  k. 

Now  measure  b  from  the  line  of  nodes. 

b  =  cos(d>  -  fl)  cos  Si'  -f  sin(<£  -  Cl)  cos  <5j '  -f  sin  <5k. 
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a  and  c  are  therefore 


a  =  —  sin  *j'  +  cos  ik 
c  =  cosij'  +  Bintk. 

Portion  of  the  Reference  P lone 


Since  a  and  b  are  perpendicular,  a  ■  b  =  0,  therefore 

tan  <5  =  tan  i  sin («^>  -  fi). 


Continuing 


results  in 


b  •  c  =  |bj  |c|  cos(0  -  rr/2) 

sin (<^  -  H)  cos  6  cos  i  -f  sin  6  sin  t  =  sin  0. 


(5.1) 


EE 


And  using  the  relation  from  equation  (5.1)  results  in 

sin£  =  sin#  sint.  (5.2) 

Finally,  from  c  x  b,  the  following  relationship  is  obtained 

cos  6  =  cos  0 sec  (0  —  H).  (5.3) 

(5.1),  (5.2),  and  (5.3)  are  the  required  angle  relations. 

D.  THE  POTENTIAL 

In  Chapter  III,  a  simplified  gravitational  potential  GM/r  was  used  with  the 
assumption  that  the  attracting  body  was  spherically  symmetric.  This  simplified 
potential  caused  the  satellite  to  move  in  conic  orbits.  As  has  been  stated,  the 
planets  are  not  spherically  symmetric  but  are  bulged  at  the  equator,  flattened  at 
the  poles,  and  are  generally  asymmetric.  An  expanded  expression  for  the  gravi¬ 
tational  potential  will  be  developed  in  this  section.  The  final  expression  for  the 
potential  is  subject  to  the  assumption  made  in  Section  B  of  this  chapter  that  the 
planet  about  which  the  satellite  revolves  is  approximately  an  oblate  spheroid. 

Now,  regardless  of  the  nature  of  an  attracting  body,  the  potential  V  must 
satisfy  one  of  the  following  differential  equations.  For  regions  within  the  attracting 
matter 

VJP  =  4npG  where  p  is  the  density  (5.4) 

which  is  Poisson’s  equation. 

For  regions  outside  attracting  matter 

V2V  =  0  (5.5) 


which  is  Laplace’s  equation. 


This  discussion  will  be  restricted  to  Laplace’s  equation  which,  if  V  is  written 


as  a  function  of  the  spherical  coordinates  (r,fl,  <f>),  may  be  written  as 


}_d_(  ,dV\  1  d  (  dVA  1  d2V 

r2  dr  ^  dr  )  +  r2  cos  6  d6  \  d6  J  r!  cos2  6  d<j> 2 

Any  conservative  force  field  F  can  be  written  as  the  gradient  of  a  potential. 

Equations  (5.5)  and  (5.6)  state  that  div  (grad  V)  is  zero.  The  assumption  that  F 

is  rotationally  symmetric  means  that  dV/d<t>  is  zero.  Therefore  a  solution  of  (5.6) 

is  sought  that  is  a  product  of  a  function  of  r  alone  and  a  function  of  6  alone: 


V(r,6)  =  h(r)  ■  f2(S). 


By  multiplying  (5.6)  by  r2,  dividing  by  /i(r)/2(6),  and  rearranging  terms,  the 
following  result  is  obtained. 


1  d 

T\dr  '  r" 


dr 


cos  6 


df± 

dd 


/2  cos  6  86 

Since  the  left  side  is  a  function  of  r  alone  and  the  right  side  is  a  function  of  6 
alone,  both  members  are  equal  to  a  constant,  say  n(n  +  1).  f\  must  satisfy 


d_ 

dr 


dr 


=  n(n  +  l)/i 


which  has  the  solution 


fu  =  Ar~n~l  -f  Brn. 


It  is  desired  that  this  function  be  zero  at  r  =  oo,  and  that  it  be  analytic  and 
single-valued  there.  Therefore  n  is  chosen  to  be  an  integer  greater  than  zero,  and 
B  most  equal  zero.  Therefore 

/i  =  Ar""-1. 


The  equation  for  /2  is 


d_ 

dS 


+  n(n  +  l)  cos  6/2=0 
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which  may  be  rewritten  as 


dsm6 


Jc> 


(1  —  sin  S)—j 


+  n(n  +  l)/2  =  0. 


(5.7) 


dsin£j 

This  is  known  as  Legendre’s  equation.  Using  the  method  of  Frobenius,  a  series 
solution  to  Legrendre’s  equation  is  found  to  be 

ft  =  CPn(sin£). 

Therefore  the  solution  to  (5.6)  may  finally  be  written  as 

V(r,5)  =  .Ar~n-1Pn(sin6) 
where  C  is  chosen  as  1. 

There  are  many  ways  to  consider  the  polynomials  Pn.  Rodrique’s  formula  for 
P„(sin  6)  is 


so  that 


P„(sin  6)  =  — -jj  f*--  -  (sin2  6  -  l)n. 
v  ‘  2nn!  d(sm  <$)n  v  ; 


P0  =  1,  Pi=sin6,  P2  =  ^(3sin2^  —  1),  etc. 


The  general  solution  of  (5.6)  is  then  written  as 

I'M)  = 


(5.8) 


n=0 


The  above  mathematical  derivation  of  an  expression  for  the  potential  may 
be  enhanced  by  a  more  physical  formulation.  Referring  to  Figure  5.5,  let  a  unit 
mass  m  be  placed  at  a  point  P  which  is  a  distance  r  from  the  center  of  mass  of 
a  bounded  distribution  of  total  mass  M.  Let  dm  be  an  element  of  the  mass  at  a 
distance  £  from  O.  Then  the  potential  at  P  due  to  dm  is 


dV 


—  K2dm 

(£2  -r  r7  -  2r£cos0)1/2 
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and  the  total  potential  is 


Figure  5.5:  Potential  at  an  Exterior  Point  Due  to  an  Irregular  Mass 


The  denominator  of  the  right  side  of  (5.9)  may  be  expanded  such  that 
(l2  +  r2  -  2rlcos0)-1/2  =  ^,(sin6) 

where  S  =  tt/2  -  6. 

Therefore  the  following  may  be  written 


V  =  f  dm-  ~r  (  £  sin  6dm  +  ~  f  l2  dm 

r  Jm  r2  Jsf  2 rs  Jm 

—-—r  f  l2  sin2  6  dm  +  . . . 

2  r3  Jm 

A  description  of  the  physical  significance  of  each  of  the  above  integrals  can  be 
made.  The  first  integral  is  the  total  mass,  the  second  is  the  first  moment  about 
an  axis  through  O  perpendicular  to  OP  and  is  zero  when  the  origin  is  chosen  as 
the  center  of  mass.  The  third  integral  is  the  moment  of  inertia  about  the  origin, 
and  the  last  integral  may  be  written  as 


-  t2  cos2  6)dm  =  I0  -  I 


where  /o  is  the  moment  of  inertia  about  the  origin  and  I  is  the  moment  of  inertia 
about  the  line  OP. 
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The  potential  is  then 


V 


K*M 


^(2/o-37). 


Now  the  first  term  is  the  potential  due  to  a  homogeneous  solid  sphere.  The  second 
term  arises  from  the  departure  of  the  mass  M  from  spherical  shape. 

With  the  physical  description  established,  (5.8)  may  be  written  as 


V  =  - 


GM 


(5.10) 


Note  that  only  Po  and  P2  remain. 

An  advantage  of  (5.10)  is  that  it  can  immediately  be  written  down  once  axial 
symmetry  has  been  assumed,  and  any  suitable  experiment  (such  as  an  orbit  of  a 
satellite  about  the  Earth)  can  be  used  to  determine  J2  (or  J3,  J4,  etc.;  for  general 
potential). 

E.  THE  EQUATIONS  OF  MOTION 

Referring  again  to  Figure  5.1,  the  components  of  the  velocity  of  a  satellite  in 
the  coordinate  system  are: 

V  =  VT  +  +  Vi 


or 


dr  d<t>  .  d6 
v  =  —  +  r  — -  cos  Hr  —  . 
dt  dt  dt 


An  expression  for  the  kinetic  energy  ( T )  is 


(drY  2  2  c  (d4>Y 

2T  =  —  +  rM  —  +  r2  cos2  6  [-f-  . 

\dt  \  dt  )  \dt 


(5.11) 


The  equations  of  motion  may  be  written  using  Lagrange’s  equation 
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(5.12) 


d  d(T  -  V)  d(T  -  V) 


=  0  q  =  r,8,  or  <j> 


dt  dq  dq 

where  q  =  T  is  the  kinetic  energy  and  V  is  the  potential  energy. 

Applying  (5.12)  to  the  expression  for  T  (5.11)  and  V  (5.10)  results  in 


dV 

dr 


dPr  ( d6\ 2  ic(d<t>V 

—  —  r  I  —  )  —  r  cos  6  1  —  1  =  — 

dt2  \dt )  \  dt ) 

d  (  j dS\  2  .  .  c  (diV  dV 

dt  \  dt  J  \dt  J  dS 

—  ( r 2  cos2  6^j-  \  =  0 

dt  \  dt  J 


(5.13) 

(5.14) 

(5.15) 


The  last  of  these  equations  gives  an  integral  which  can  be  used  to  eliminate 
time  from  the  system,  i.e., 

.  d(f> 


Let  r  =  p/u  then 


r2cos2  8—  —  h  cos  i0. 
dt 


d  _  h  cos  t0 u2  sec1 6  d 
dt 


P 2 


d<j> 


(5.16) 


Starting  with  equation  (5.14)  and  applying  (5.16)  to  the  first  term  results  in 


dt  \  dt  j  dt 


d  I  r2h  cos2  i0scc26\  d6 
p2  )  d<f>' 


and  continuing  to  apply  (5.16)  so  that  (5.14)  finally  becomes 


d 2 1 '  ri  8  dVr2cos28 

-1-  tan  8  =  —  • 


d4>2 


d8  h2 


COS'  to 


or 


d2  tan  8  r  sin  8  cos3  8  .  _  , 

+  tan  8  -  - — - (2Ju) 


d4>2 


COS'  »o 


(5.17) 
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where 

'  ■  i  (*a 

p  =  h2/GM. 

It  is  desired  to  get  the  left  side  of  (5.17)  in  terms  of  the  independent  variable 
6.  To  begin,  the  right  side  of  the  angle  relationship  (5.1), 


tan  6  =  tan  t  sm(4>  —  fl), 


is  substituted  into  (5.17). 
Then  using  (5.3) 

(5.17)  becomes 


cos  8  =  cos  6  sec(4>  -  fl), 


/  W\  2  sec2 1  cos  Of  t'  tan  » cos  6  l  Cl ’<j>"  -  Cl" 

l  4>' )  cos  8  l 1  <j>' )  <j>'  cos  6  l  4>n 


-2  sec2 1  tan  8 


<t>'  I  sin 


an  8  (  i"  - 

i  cos  i  \  <))'2 


sin  8  cos3  8 
cos2  to 


(2  Ju), 


(5.18) 


where  primes  denote  differentiation  with  respect  to  6. 


Later,  4>  will  be  eliminated  from  (5.18),  but  first  the  remaining  equation  of 
motion  (5.13)  will  be  rewritten  in  a  form  like  (5.18). 

To  rewrite  equation  (5.13)  in  a  useful  form,  the  independent  variable  is  again 
changed  by  use  of  (5.16).  Multiplying  the  expression  by  (d>')2  and  noting  that 

(tan2  <5)'  =  2  tan  8  sec2<5  — 
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results  in: 


u 


+  u 


,,  2  (tan£)'2 

‘  C0S  ‘  +  M> 


+ 


u'(tan2  5)' 
sec2  6 


<t>' 


_  <p'  cos4  6  dV 
h2tt2  dr 

Evaluating  ^  and  multiplying  (5.19)  by  results  in: 


(5.19) 


u"  +  u  = 


cos2  i 


COS*  *0 


jr-[l  +  J  u2(l  —  3  sin2  £)]  +  u" 


1  - 


(l  +  n'~4  —  *'  cos  6  sin  6  tan  t 

\  cot  I 


u'  cos2  x  sec4  6 


All 

2  sin  6  cos  6 (tan  6)' - 

<P 


—  u 


cos2 1 


sec2  6  + 


(tan  <5)'2  cos2 
ft* 


-  1 


(5 


Equations  (5.18)  and  (5.20)  are  the  equations  of  motion.  Before  proceeding 
to  solve  these  equations,  it  is  necessary  to  remove  <f>.  Combining  equations  (5.1)  - 
(5.3)  results  in  the  following  expression: 


tan(d>  —  H)  =  cos  i  tan  0. 

Differentiating  this  expression  with  respect  to  0  results  in 
,,  cos  i  .  *'  sin  6  cos  6  sin  t 

<i>  —  - +  n  —  - — - 

cos*  c  cos*  o 


(5.21) 


or 


1  cos2  6  1 

<t>'  cos  i  1  +  f] sin  $  cos  $  tan  t 

COB  1 

From  (5.2),  sin2  6  —  (sin  6  sin  :)2  or 


(5.22) 


cos2  6  —  1  -  (sin  6  sin  x)2. 
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Substitution  of  this  expression  plus  (5.21)  and  (5.22)  in  (5.18)  and  (5.20) 
results  in  completely  general  equations  of  motion  in  terms  of  the  dependent  vari¬ 
ables:  t,  fl,y;  and  the  independent  variable  0. 


6  Sil1  *  (S)lsinS*  3  Sin  dfi  COS  6  8in2  6 
cos  t  cos  t 

3  sin2  * jj  cos2  0  sin  6  ^  6  sin3 1  (^)  *  cos2  6  sin  0 


cos  t 


cos  t 


+ 


3  sin2 t 4jl  cos2  0  sin  6  2  sin2  sin  6 


cos  » 


cos  t 


sin  0 


.  (ft  d  n  .  ^ ..  1.3 

+  3  cos  t—  — —  sin  6 - - 

dO  d6 2  cos  t 

3sin3i(f)2sinfl  .  .  .  ( dU\2  .  „ 

+  - 7 - (-  6  cos  t  sin  1  —  sin  0 

cos  1  J 

3sint(f)2sinJ  ^Sin2,gif  sinfi 


cos  t 


ld6i  d$ 

cos  t 


,d2i  dn  ,  jkdt  sin  0 

+  3  cos  t  — -  —  sin  6  -  &  -de  -- — 
dO 2  dO  cos  » 


dn 

—  2  sin  t—  sin  6 
dO 


+ 


d2i  . 

5«JS'nS-  cos, 
12sin*igf  cos»« 
cos  l 


6  sin  »^2  sin  6  3  sin3 cos3  0 


de  de* 
cos  l 


.  .dn  d2n 


—  3  cos  i  sin  i—— 


dO  dO2 


cos  6 


.d2 Q  2  sin2 1 ^  cos  6 

+  sin  t-—  cos  6  + - - 

do2  cos  1 

.di  dn  „  cos  0  di 

+  12  cos  1  —  —  cos  0 - . - 2—  cos  0 

dd  dO  cos  t  dO 

2  cos3  i  sin  iJu  sin  0 

COS2  tr, 


(5.23) 


and  equation  (5.20)  is: 
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+ 


+ 


+ 


+ 


+ 


sin 


3sin4»(f)*ucos4fl 
cos2  i 

6  sin  i£4g%cos20 


cos2  t 


■_3  ■  di  x in  du  /] 

Sln  *d6  dt  dt  COS  _ 

cos2  i 


<t*  dt  dt 

cos2  t 


.  .dt' dfi  du  i/t  3  sin  cos2 6 

3  sin  t—  —  —  cos1 0 - ^ - 

d6  d0  d6  cos  t 


-  3  sin=  ii^X u  cos1  S  +  i"’1' ife C°S’  *  -  iitpl 
\d0  J  cos  t  cos2 t 


ydn&u  odll  d'il  du 

*  J/12  “  Jit  J«2 


'  de  dp 
COS  1 


in  i^n  du 


dt  dt*  de 
cos2 t 


+ 


xpn  du 

dt*  dt 


COS  t 

2  sin  i£e§fe  .  <*fi  du  sin  t£^ 

- de  de  de  -  sin  t  —  —  —  + - 

dfl  dd  dO 


cos2  * 


cos  t 


cos 


2  iJ  u2 


COS2  t’o 


+ 


cos2  t 
cos2  to 


(5.24) 


Equations  (5.23)  and  (5.24)  can  be  greatly  simplified  if  only  a  first  order 
approximation  to  the  equations  of  motion  is  desired.  The  approximate  equations 
to  (5.23)  and  (5.24)  are  respectively: 


.  .dH  . 

2  sin  t  —  sin  6 
dO 


d2t'  .  „  .  ,d2n  „  dt 

—i  sin  0  +  sin  t— —  cos  6  -  2—  cos  6 
ddL  dd 1  dd 

2  cos3 t  sin  t J  u  sin  6 
cos2 10 


(5.25) 


—  +  u=  sin  l^sin(2g) 
d02  cos  t 

2  sin2l’f  £rsin2e  _  sin2tg£^sin2fl  _  sin  t'^sin2fl 

cos  t  cos2 t  cos  1 

3  cos2 t  sin2 1  Ju 2  sin2  0  2  sin2 t  ^  cos  0  sin  0 

cos2 t  cos  t0 

sin  t  ^  ^  cos  0  sin  0  2  sin  t^u  cos  0  sin  0 

cos  t  cos  i 
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(5.26) 


3  sin  ijjjjcos2  6  2  sin2  i^u  cos2  0 


cos  t 


cos  t 


+ 


+ 


nd£i£u  d?n  du  •  r-_J  •  7  ,.2 

~  d»  M*  _|_  d$2  dS  _j_  Sln  1  dS  dt  +  COS  tJU 

cos  i  COS  t  COS  t  COS2t'o 


cos2 »' 


COS'*  *0 

Equations  (5.25)  and  (5.26)  will  be  used  in  the  initial  analysis;  the  general 
equations  will  be  used  for  second  order  calculations. 


F.  THE  INITIAL  CONDITIONS 

It  is  desired  that  the  solution  to  the  equations  derived  in  the  last  section 
be  an  accurate,  long  term  predictor  of  the  satellite’s  motion.  In  fact,  as  long  as 
the  oblateness  perturbation  remains  the  dominant  disturbing  force,  the  solution 
should  be  valid  for  close  to  1000  revolutions.  However,  before  the  solution  can  be 
used  for  prediction,  a  set  of  initial  conditions  must  be  determined. 

The  subject  of  orbit  determination  represents  a  separate  discipline  within 
celestial  mechanics,  and  a  discussion  of  the  various  techniques  used  to  determine 
an  orbit  is  outside  the  scope  of  this  analysis.  It  is  sufficient  therefore,  to  state  that 
the  purpose  of  orbit  determination  is  to  find  the  orbital  elements  of  a  satellite  from 
reduced  observational  data.  A  set  of  observations  will  determine  the  osculating 
elements  at  time  t0.  As  was  noted  in  Chapter  4,  these  elements  will  change,  and 
at  a  new  set  of  osculating  elements  may  be  calculated.  For  the  purpose  of  this 
analysis  any  observed  set  may  be  designated  as  the  osculating  elements  at  t0  and 
thus  prescribe  the  initial  conditions.  Staged  mathematically,  the  initial  conditions 
are  at  t  =  t0: 


a(l  -  e2)  dr  _  o(l  -  e2)e  sin(#o  —  ^’) 

1  +  e  cos(0c  -  u.’)  d8  (1  +  f  cos(6q  -  w)2 
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where  a  =  oq,  e  —  to,  u>  =  u>o, 

and  t  =  to  fi  =  H0  y  =  0o-u>  (5.27) 
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VI.  THE  PERTURBATION  PROCEDURE 


A.  PRELIMINARIES 

The  perturbation  method  used  in  this  analysis  is  a  variation  of  tne  technique 
known  as  the  Method  of  Strained  Coordinates.  The  motivation  for  the  use  of  this 
technique  is  the  subject  of  this  section. 

In  perturbation  theory,  the  quantities  to  be  expanded  can  be  functions  of  one 
or  more  variables  besides  the  perturbation  coordinate.  An  asymptotic  expansion 
of  f(6\J )  in  terms  of  the  asymptotic  sequence  em(J)  is, 

f{B\J)  ~  £  am(0)em(J)  as  J  — ►  0  (6.1) 

m=0 

where  6  is  a  scalar  (or  vector)  variable  independent  of  J  and  the  coefficients  am 
are  functions  of  6  only.  This  expansion  is  said  to  be  uniformly  valid  if, 

=  £>  m(«)emU)  +  *n(M) 

m= 0 

Rn(6;J)  =  0(£„(J)). 

For  these  uniformity  conditions  to  hold,  am(6)cm(J )  must  be  small  compared  to  the 
preceding  term  for  each  m.  Each  term  must  be  a  small  correction 

to  the  preceding  term  regardless  of  the  value  of  6. 

Unfortunately,  it  is  the  rule  rather  than  the  exception  that  expansions  like 
(6.1)  are  non-uniformly  valid  and  break  down  in  certain  cases.  A  case  of  particular 
interest  is  the  presence  of  secular  terms  such  as  6n  cos  6  and  6n  sin  6  which  make 
fm{8)/ fm-i  (0)  unbounded  as  6  approaches  infinity. 

To  illustrate  how  secular  terms  arise  in  the  solutions  to  differential  equations, 
reference  is  made  to  one  of  the  equations  of  motion  derived  in  the  preceding 
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chapter.  Equation  (5.26)  has  the  form, 


u"  +  u  =  i  +  j/(u,u\u'\  n',n") . 

When  expanded,  the  right  side  will  be  linear  combinations  of  trigonometric  terms 
and  constants.  The  presence  of  terms  of  the  form  cos(0),  sin(0),  cos(3y  —  26), 
sin(2y  —  6),  etc.,  on  the  right  side  will  produce  secular  terms.  For  example,  the 
second  order  differential  equation, 

u"  +  u  =  cos  6 


has  the  solution, 


6  sin  6  cos  6  . 

u  =  - — — — h  Ki  sin  6  -t-  Ki  cos  6. 


Note  that  the  presence  of  the  term  6  sin  6  will  result  in  unbounded  solutions  for 
u  as  6  approaches  infinity.  Since  in  this  representation  u  is  the  reciprocal  of 
the  distance  from  the  center  of  the  planet  to  the  satellite,  the  physical  effect  of 
the  secular  terms  would  be  to  produce  in  r  periodic  terms  with  large  amplitude 
variations,  a  situation  that  certainly  does  not  correspond  to  physical  reality. 

A  technique  for  dealing  with  terms  that  produce  secular  terms  is  to  eliminate 
them  from  the  right  side  of  the  differential  equation.  The  Method  of  Strained 
Coordinates  is  a  perturbation  technique  designed  for  removing  secular  terms.  To 
illustrate,  it  is  recalled  that  the  solution  to  the  two  body  problem  is 


l/r  =  u  =  1  +  e  cos (6  -  u>),  where  p  =  1. 

In  the  Method  of  Strained  Coordinates,  the  6  coordinate  is  strained  by  introducing 
a  new  variable  y  —  /  (0)  =  6  —  u>  -f  Jy\6  +  ...  As  will  be  shown  in  the  next  section, 
choosing  the  correct  value  for  y j  will  insure  secular  terms  do  not  arise  in  the  first 
order  solution  for  u. 
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While  the  above  method  removes  secular  terms  to  first  order,  it  is  not  ade¬ 
quate  for  dealing  with  secular  terms  that  arise  to  second  order  in  the  equations  of 
motion.  That  it  is  important  to  remove  second  order  secular  terms  is  shown  by 
the  following  equation: 

f(r)  =  1  +  t  cos(0  —  w)  -1-  J(cos  8  +  . . .) 

+  J2(0  cos (6  —  tjj)  +  . . .)  +  J3(02  cos(0  -  w)  +  ...)  +  .. . 

Note  that,  although  the  terms  through  order  J  are  bounded,  as  6  — +  oo,  the  J 2 
and  J3  terms  grow  without  bound  and  dominate  the  right  side.  In  this  present 
analysis,  6  has  an  upper  bound  of  (27r)l03;  however,  the  effect  of  secular  terms 
remains  since  J26,JS$2,  etc.,  are  all  of  order  J  as  6  — >  (27r)10s.  An  infinite  series 
would  have  to  be  retained. 

In  this  present  work,  three  additional  techniques  had  to  be  devised  to  deal 
with  secular  terms  to  second  order.  These  techniques  will  be  discussed  in  the  next 
section.  The  perturbation  method  therefore  is  not  strictly  the  Method  of  Strained 
Coordinates,  but  a  variation  of  that  method. 

The  basic  steps  in  the  perturbation  procedure  are  as  follows: 

1.  The  dependent  variables  and  independent  variables  are  expanded  as  func¬ 
tions  of  a  small  parameter  ( J ). 

2.  The  variables  are  then  substituted  into  the  equations  of  motion,  and  the 
equations  are  solved  consecutively.  Each  solution  yields  a  more  exact  ex¬ 
pression  for  the  appropriate  variable. 

The  process  is  carried  out  through  second  order  to  demonstrate  that  all  secular 
terms  may  be  eliminated  and  that  the  solutions  are  bounded.  The  following  sec¬ 
tion  highlights  the  calculations  involved  in  the  process  and  shows  the  first  order 
equations  and  solutions.  The  second  order  expressions  are  long,  and  their  display 


in  this  context  would  not  contribute  to  tne  analysis  since  only  a  few  specific  terms 
are  relevant.  The  complete  expressions  are  contained  in  Appendix  B;  they  will  be 
referred  to  during  the  course  of  the  analysis. 


C.  SOLVING  THE  EQUATIONS 

1.  First  Order  Approximation  for  i  and  H 

The  equations  to  be  solved  are  (5.25)  and  (5.26).  Since  the  right  side 
of  these  equations  are  analytic  functions  of  J,  it  is  reasonable  to  expect  that  the 
solution  to  u  will  be  arbitrarily  close  to  the  two-body  solution,  1  +  e  cos(0  —  w), 
when  J  is  sufficiently  close  to  0.  Likewise  t,  fl,  and  y  will  be  arbitrarily  close  to 
j'o,  Do,  and  0  —  u,  respectively,  when  J  is  close  to  0.  This  assumption  amounts  to 
letting 


1  +  e  cos  y  +  Jm  +  J2u2  +  . . . 

(5.2) 

6  —  u)  -f~  Jy\  J2y2  +  . . . 

(5.3) 

t'o  -f  Jtj  +  J^i2  +  . . . 

(5.-:) 

fio  +  Jfii  +  «/2fi2  +  •  •  • 

(5.5) 

The  first  step  in  the  solution  of  the  equations  of  motion  is  to  substitute 
(6.2)  -  (6.5)  into  equation  (5.25)  and  equate  the  coefficients  of  J .  The  result  is 

(—2n[  sin  6  +  f2"  cos  6)  sin  i0  —  2t’x  cos  6  -  i"  sin  6 
=  2  cos  »0sin  t0sin  0(1  +  e  cos  y)  .  (5.6) 


Sin  t  and  cos  t  have  been  replaced  in  the  above  equations  by  their  approximations: 
sin  t'o  and  cos  i0.  These  are  valid  approximations  since  H",  O',  t",  and  are  all  of 
order  J. 
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Equation  (6.6)  is  a  linear  differential  equation  with  two  unknown  func¬ 
tions.  Terms  of  the  form  cos  0,  sin  6 ,  etc.,  on  the  right  side  of  this  equation,  and 
the  more  general  equation  (5.23),  will  produce  secular  terms  in  t  and  fl.  It  should 
be  recalled  that  these  same  terms  cause  secular  terms  in  u.  Note  that  (6.6)  con¬ 
tains  a  sin  6  terms  on  the  right  side.  There  would  be  a  need  to  eliminate  this 
term  were  it  not  for  the  conditions  placed  on  H  and  t  by  the  definition  of  the 
reference  plane.  Specifically,  fi,  which  governs  the  rotation  of  the  reference  plane, 
must  be  expressible  as  an  unbounded  (secular)  term  plus  bounded  periodic  terms, 
and  *  must  be  bounded.  The  fact  that  Cl  must  contain  a  secular  term  negates  the 
requirement  to  eliminate  the  sin  6  term.  Therefore,  a  solution  of  the  form 

fti  =  a\6  4-  a.2  sin  y 

*i  =  Pi  cos  y  (5.7) 

is  assumed. 

By  substituting  (6.7)  into  (6.6)  and  equating  coefficients,  the  following 
solution  for  i  and  fl  is  obtained: 

t"i  =  — 2/3  e  cos  t'o  sin  t’c  cos  y 

Cli  =  —6  cos  t'o  —  4/3  e  cos  i’osin  y  (5-8) 

2.  Second  Order  Secular  Terms  in  i  and  fl 

Equation  (6.8)  satisfies  (6.6)  to  order  J .  The  next  step  in  the  process 
is  to  substitute  (6.2)  -  (6.5)  and  (6.8)  into  (5.26)  and  solve  for  uj.  However, 
proceeding  with  (6.8)  in  its  present  form  will  lead  to  secular  terms  in  second 
order.  A  brief  paragraph  will  explain  why  secular  terms  are  anticipated. 

Success  in  using  perturbation  methods  requires  a  certain  a  priori  knowl¬ 
edge  about  the  nature  of  the  particular  problem  one  is  trying  to  solve.  There  is  a 
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certain  amount  of  trial  and  error  involved  in  the  process.  For  instance,  the  success 
of  the  Method  of  Strained  Coordinates  depends  on  the  knowledge  that  a  secular 
term  will  arise  in  the  first  order  solution  to  (5.26)  and  that  a  mechanism  (the 
strained  coordinate)  can  be  put  in  place  one  step  ahead  in  order  to  eliminate  the 
secular  term.  The  second  order  secular  terms  were  found  by  this  same  trial  and 
error  process. 

Returning  now  to  equation  (6.8).  it  was  discovered  in  this  analysis  that 
the  problem  terms 


(-45ss  +  28s) 
24 


J2e2c  sin(2y  —  6) 


(6.9) 


and 


^  J2sse2c  sin(2y  —  30) 
8 


(where  s  =  sin  i0,  c  =  cos  t0 ) . 

appear  on  the  right  side  of  (5.23)  when  it  is  expressed  to  order  J 2.  The  appearance 
of  these  terms  that  give  rise  to  secular  terms  was  unexpected.  They  were  not 
reported  by  Brenner  because  the  authors  assumed  a  small  eccentricity  and  dropped 
J2e2  and  higher  terms. 

There  were  several  failed  attempts  in  dealing  with  these  new  terms  be¬ 
fore  an  effective  measure  was  discovered.  First,  an  investigation  was  made  into 
the  effect  of  retaining  the  secular  terms. 

The  secular  terms  produced  by  (6.9)  are: 

sin(2y  -  2g) 

Q2  =  ^  ^ — 'c:c  6  cos(2y  —  2d) . 


As  was  demonstrated  earlier,  when  6  — ►  (27r)lOs  these  terms  are  of  order  J  and 
must  be  retained  in  the  first  order  solution.  As  the  solution  progresses,  these  terms 
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will  continue  to  produce  secular  terms  with  coefficients  J563,  J*6S,  . . Jn0n~1.  A 
convergent  series  representation  for  these  terms  was  not  found. 

The  remaining  alternative  was  to  alter  the  form  of  (6.7)  to  eliminate 
secular  terms.  If  terms  with  arbitrary  coefficients  can  be  found  which  when  added 
to  (6.7)  produce  terms  of  identical  harmonics  to  the  same  order  as  (6.9),  the  co¬ 
efficients  may  be  chosen  so  that  all  terms  with  these  harmonics  are  eliminated. 
Essentially,  the  only  terms  that  may  be  added  to  (6.7)  are  terms  which  satisfy  the 
original  differential  equation,  equation  (5.23).  Therefore,  one  may  add  homoge¬ 
neous  solutions  of  the  differential  equation  or  arbitrary  constants. 

Adding  an  arbitrary  constant  to  (6.7)  has  the  following  effect  on  higher 
order  terms.  From  equations  (5.23)  and  (5.24),  it  would  appear  that  only  the 
derivatives  of  i  and  ft  enter  into  the  higher  order  calculations,  and  therefore  arbi¬ 
trary  constants  would  be  eliminated.  This  is  true  for  ft,  but  not  for  i.  To  explain, 
it  is  recalled  that  the  approximations  sin  i  —  sin  t0  and  cos  i  -  cos  i0  were  valid 
for  the  first  order  approximation  of  tj  and  ftj.  This  approximation  was  valid  since 
all  terms  were  of  order  J.  However,  in  the  calculation  for  U!  there  is  a  term 
cos2  ij  cos2  t‘o  that  is  of  order  1  (Eq.  (5.26)),  therefore  a  better  approximation  for 
:  is  required.  Using  (6.7)  and  the  Taylor  series  expansion  for  cos  i  (keeping  only 
terms  of  order  J)  results  in 

cos  t  =  cos(i0  -f  /(J))  =  cos  i0  -  sin  t0 f(J) 

=  cos  t’o  —  sin  *o(*i) 

Adding  an  arbitrary  constant  Ka  to  (6.8)  results  in 

cos  i  =  cos  i'o  4-  sin  t'o(2/3  J  e  cos  t0sin  t0cos  y  +  Ka) 

(A  similar  expression  is  required  for  sin  i) 
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Therefore,  a  constant  added  to  (6.8)  will  alter  the  form  of  subsequent 
terms,  but  not  in  the  form  required  to  eliminate  secular  terms  in  t2  and  fl2.  It  will 
be  shown  that  a  constant  like  Ka  will  produce  a  term  of  the  form  KaJ  sin2  »0  in  ui 
and  terms  of  the  form  KaJ2c  sin2 10  cos  i0  cos  y  in  »2,  and  KaJ2e  sin2  t0  cos  t0  sin  y 
in  n2.  Although  they  cannot  be  used  to  eliminate  secular  terms,  the  constants 
will  be  needed  to  satisfy  the  initial  condition,  therefore  it  is  essential  that  they  do 
not  produce  irremovable  secular  terms  to  higher  order. 

The  next  alternative  for  the  elimination  of  secular  terms  is  to  add  a 
solution  to  the  homogeneous  equation  of  (5.23).  While  it  will  be  shown  that  this 
technique  is  successful  in  eliminating  secular  terms  in  u2,  it  did  not  succeed  in 
removing  them  in  t2  and  fl2.  It  failed  because  the  homogeneous  terms  produced 
new  secular  terms  to  higher  order.  Many  various  combinations  of  homogeneous 
solutions  were  added  to  (6.7),  but  all  attempts  along  this  line  only  complicated 
the  problem. 

An  answer  to  the  question  of  how  to  eliminate  the  secular  terms  was 
suggested  in  a  report  by  Weisfield  [Ref.  31]  on  polar  orbits.  When  faced  with  the 
problem  of  eliminating  secular  terms  from  his  equation  for  A2,  Weisfield  added  a 
term  like  cos(2y  -  26).  Now  to  the  particular  order  to  which  one  is  working  this 
term  acts  like  a  constant  in  the  derivative,  e.g.,  let, 


y  —  6  +  JO 


then 


—  (cos(2y  -  26))  =  sin(2y  -  20){2  -  2(1  +  J))  =  O(J). 
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The  attempt  was  then  made  to  apply  Weisfield’s  technique  to  this  more 
general  problem.  (6.7)  would  then  become 

j'i  =  -2/3e  cos  t0sin  t0cos  y  +  Kx  cos(2y  —  26)  +  K2 

fii  =  6  cos  t'o  -  4/3ecos  i0sin  y  +  jfjsin^y  —  26)  +  K4  .  (6.10) 

where:  Kx  =  -Kxe7cs  K2  =  Kx  cos(2o>)  +  K2 
Ki  =  -Kse7c  K4  =  -Ks  sin(2w)  +  K4 

(It  is  noted  that  (6.10)  still  satisfies  (6.6)  to  order  J) 

The  added  harmonic  terms  produce  identical  harmonics  as  (6.9)  to  sec¬ 
ond  order  and  no  other  new  secular  terms.  Thus,  they  allow  for  the  elimination 
of  these  problem  terms. 

3.  First  Order  Approximation  for  u 

With  valid  expressions  for  t’i  and  Clx  established,  the  procedure  may  now 
be  continued  with  the  calculation  of  ux.  Substitution  of  (6.2)  -  (6.5)  and  (6.10) 
into  (5.26)  and  letting  sin  *0  =  s  and  cos  »o  =  c  results  in  the  following  equation 

d2ux  5e2s2  cos(2y  -f  26) 

- +  ux  =  2e  cos  y  yx  +  -- 


d6 2  1  . * 8 


+  2e2Kxs2  cos(2y  —  26)  + 


e2s2  cos(2y  —  26) 


8 


5es2  cos(y  26)  lle2s2  cos(2y)  5e2cos(2y) 


2  3e2s2cos(20)  s2 cos (26) 

-  bes  cos  y  -f  4e  cos  y  -i - 1 - 

4  2 

22  2  17e2s2  5s2  7c2 

—  2e  Kxs  cos(u.')  —  2K2s - — - -  I — —  1-  1 

1 1  l  o 


(6.11) 


In  the  above  equation,  the  cos  y  terms  will  produce  secular  terms  in  ux.  A  choice 
of  yi  =  (5s2  -  4)/2  will  eliminate  that  possibility,  y  becomes 


y  =  6  -  <*>  *  J  6  i^-s2 


J2y:6. 
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Equation  (6.11)  becomes 


.  rj  5eV(2y  +  20)  a  *  «flv 

+  u i  — - - - t-2eiiis  cos(2y  —  20) 

e2s2  cos(2y  —  20)  5es2  cos(y  +  20)  lle2s2  cos(2y) 

+  8  +  3  4 

5e2cos(2y)  3c2s2cos(20)  s2cos(20) 

+  -  +  -  +  2 

*i«2  7*2 

—  2e2ifi«2  cos(2w)  -  2 K^s1 - 1 - 1-  1.  (6.12) 

12  2  6 

To  solve  (6.12),  a  solution  of  the  form 

Ui  =  ai  cos(2y  +  20)  +  c*2  cos(2y  —  20)  +  0:3  cos(y  +  20)  +  a4  cos(2y) 


a5  cos(20)  +  a0 


(6.13) 


is  assumed. 

Equation  (6.13)  is  substituted  into  (6.12)  and  the  coefficients  of  like  harmonics  are 
equated.  The  coefficients  are: 
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Equating  (6.13)  for  results  in: 


4.  Second  Order  Secular  Terms  in  u 

The  above  expression  for  uj  is  not  complete.  Looking  one  step  ahead 
it  is  known  (by  trial  and  error)  that  the  following  problem  terms  arise  in  the 
differential  equation  for  u2. 


(Only  the  problem  terms  are  displayed.  The  complete  expression  is  equation  (B.5), 
App.  B.) 

The  harmonics  cos(3y  —  26),  cos(y  -  26),  and  cos  y  all  cause  secular 
terms  in  u2.  and  therefore  must  be  eliminated.  In  addition  there  is  another  more 
troubling  problem  with  (6.15).  That  is,  certain  terms  have  in  their  denominator 
(5sin!iV  -  4).  and  if  t’c  =  —  sin-1  '5,  then  the  denominators  are  zero.  This  in- 


clination  is  well  known  in  Astrodynamics,  and  has  been  named  the  critical  inclina¬ 
tion,  or  more  appropriately  the  critical  inclinations  since  there  axe  two:  t0  =  63°26' 
and  117°34'.  The  problem  of  the  critical  inclination  will  be  dealt  with  later  in  the 
analysis.  In  addition,  there  is  qualitative  discussion  of  the  critical  inclination  in 
Appendix  C. 

The  task  is  to  eliminate  the  three  terms  which  give  rise  to  secular  terms 
in  u2  from  the  right  side  of  (6.15).  By  inspection  it  is  seen  that  a  proper  choice  of 
j/2  will  eliminate  secular  terms  produced  by  cos  y.  It  was  then  discovered  in  this 
analysis  that  the  addition  of  a  term  of  the  form  cos (y-  26)  —  cos(y  +  2u/)  eliminated 
the  cos(y  —  26)  problem  term.  This  term  may  be  added  since  it  is  a  solution  to 
the  homogeneous  equation.  There  remained  one  term  to  eliminate;  coi>(3y  —  26). 
This  term  can  be  eliminated  by  adding  a  term  with  the  harmonic  sin(2y  —  26)  to 
yi.  In  addition,  a  constant  was  added  to  yi  to  satisfy  the  initial  conditions.  The 
complete  first  order  expression  for  y  is  now 


y  =  V  -  w  +  J  ^0  ^s2  -  2^  +  Ks  sin(2y  -  26)  +  K^j 


+  J2y76  ■ 


(6.16) 


where  K$  =  —JK$  and  J?6  =  -JKi  sin(2w)  +  Ke- 

By  adding  —  K7(cos(y  —  26)  -fcos(y  +  2u/))  to  eliminate  a  secular  term  and  K$  cos  y 

-l-A'gsin  y  to  satisfy  initial  conditions,  Uj  becomes 

e2sJ  cos(2y  +  26)  ,  T.  ,  , 

uj  = - - - -  +  2e2 K \s2  cos(2y  -  26) 

24 

eV  cos(2y  -  2J)  _  5e,’  cos(v  +  20)  _  _ 

8  24  V  ' 

,  ,  lle2s2 cos(2y)  5e2cos(2y)  .  .  . 

+  eK7  cos(y  4-  2i u)  + - — - - - -  +  K%  sin (y) 

1 Z  D 


Kg  cos(y)  - 


e2s2cos(26)  sJcos(2^)  2  2 


-  2e2K^sl  cos(2u’) 
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2ff2s2  - 


(6.17) 


17eV  5s2  7e2 

12  ~  ~2  +  6  + 

5.  Second  Order  Solution  for  i  and  Q 

With  all  terms  in  place  to  deal  with  secular  terms  the  calculations  are 
continued  by  substituting  (6.2)  -  (6.5),  (6.10),  (6.16),  and  (6.17)  into  (5.23)  to  solve 
for  fi2  and  »2,  and  to  evaluate  the  constants  Kx  and  K3.  The  result  is  contained 
in  Appendix  B,  i.e., 


(— 2ft'2  sin  6  -f  fl2  cos  0 )  sin  t0  —  2t2  cos  6  —  t2  sin  6 
=  Right  side  (B.l),  (Appendix  B). 


An  inspection  of  (B.l)  reveals  that  the  coefficients  of  sin(2y  —  6)  and  sin(2y  —  3 0) 
form  respectively  the  following  simultaneous  equations: 

(5 K3  +  10 K,  -  45)e2s3c  {4K3  +  \KX  -  7)e2s  c 

24  “  6  “  ° 


~  j  e2s3c  + 


(4Ki  -  4K3)e2sc  =  0. 


Solving  these  equations  results  in 


Kx 


K3 


15 52  -  14 
24(5sJ  -  4) 

75s4  -  120s2  +  56 
24(5s2  -  4)2 


(6.18) 

(6.19) 


Substitution  of  these  values  for  Kx  and  K3  into  (B.l)  gets  rid  of  all  sin(y  —  26) 
and  sin(2y  —  3 6)  terms.  With  an  assurance  that  no  secular  terms  will  arise,  the 
equation  for  H2  and  t2  can  be  solved. 

Before  progressing,  it  is  noted  that  (6.18)  and  (6.19)  contain  the  same 
problem  denominator  that  was  observed  in  (6.15).  In  fact,  (6.18)  and  (6.19)  are 
the  first  occurrence  of  the  critical  inclination  term  in  this  analysis.  The  term  then 
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continues  to  manifest  itself  in  all  higher  order  analyses.  The  reason  for  the  presence 
of  the  term  is  simple  to  explain.  The  constants  K\  and  Ks  were  multiplied  by  the 
derivative  of  yx  —  2))  during  the  analysis.  The  derivative  then  shows 

up  in  the  denominator  of  these  constants  when  they  are  evaluated.  It  will  be 
important  to  show  later  in  the  results  that  despite  the  occurrence  of  an  apparent 
singularity,  the  final  solution  is  uniformly  valid  for  all  inclinations,  including  the 
critical. 

With  secular  terms  removed  the  following  equation  can  be  solved  for  fl2 

and  t2. 

( — 2rij  sin  0  +  Hj  cos  6)  sin  t0  -  2i'2  cos  6  —  i2  sin  6  (6.20) 

=  Right  side  (B.2),  (Appendix  B.). 

As  was  done  in  (6.6),  a  solution  which  includes  the  harmonics  of  the  right  side  of 
(B.2)  is  assumed  for  fi2  and  t2.  This  solution  is  substituted  into  (B.2),  and  the 
coefficients  of  the  various  harmonics  are  equated.  Once  again  the  conditions  are 
that  n2  must  be  expressed  as  a  secular  term  and  bounded  periodic  terms,  while 
t2  must  be  bounded.  The  solution  is 


fi2  =  Right  side  (B.3) 

»2  =  Right  side  (B. 4).  (6-21) 

6.  Eliminating  the  Final  Secular  Terms 

To  complete  the  solution  the  remaining  constants  must  be  found.  Three 
of  the  remaining  six  constants,  y2,  K5  and  K7,  are  obtained  from  the  differential 
equation  for  u2.  It  is  recalled  that  these  constants  are  used  to  eliminate  secular 
terms  in  tx2.  There  is  no  need  to  solve  the  equation  for  u2  since  the  constants  may 
be  evaluated  from  the  right  side  of  the  differential  equation.  Therefore,  (6.2)  - 


68 


(6.5),  (6.10),  (6.16),  (6.17),  and  (6.21)  are  substituted  into  (5.24).  The  resulting 
equation  is  (B.5). 

As  in  the  previous  solution  for  O2  and  t'2,  the  secular  terms  in  u 2  arc 
eliminated  by  setting  the  coefficients  of  the  problem  terms,  i.e.,  (6.17),  equal  to 
zero  and  solving  for  the  constants.  The  results  are: 


where 


75e2s6  -  260e2s4  +  296eV  -  112e2 
192  (fs2  -  2) * 

15e2s4  —  (I4e2  -  2 )s2 
48  (fs2  —  2) 

e2  cos(2w) 

30  (fs2  -  2 


(6.22) 


y  2  = 


+ 


26es2cos(cj  -  60)  15e2s4  cos(2cj) 

3  8 

15e2s2 cos(2w)  e2cos(2w)  15e2s4  275s4  3 e2s2  61s2 

8  60  +  32  48~~  +  8  +  ~12 

I6! 

12 


With  the  above  constants  evaluated,  it  is  assured  that  the  second  order  expression 
is  free  of  secular  terms. 

7.  The  Initial  Conditions  for  y,  t,  8,  H 

The  task  now  is  to  evaluate  the  remaining  constants  by  establishing  the 
initial  conditions.  At  t  =  t0,  it  is  required  that  the  velocity  vector  of  the  satellite  in 
the  reference  plane  be  tangent  to  the  corresponding  two  body  ellipse  determined 
by  the  satellite.  In  addition  to  the  initial  conditions  established  for  r  and  ~  in 

OP 

(5.27).  it  is  recalled  that  at  t  —  to 


y  =  Or,  -  u\  i  =  t’o  and  fl  =  Ho 
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From  (6.16)  and  (6.22): 

y  =  0-u  +  j(e^s2-2^j 

(75e2s6  -  260eV  +  296e2s2  -  112e2)  .  .  .  .  .  xx\ 

-  - - - - rj - ~(sin(2y  -  26)  -  sin(2w)) 

1Q9  I  5. *2  _  91  / 


+  K6 


192  (§s2-2)‘ 


30  (§s2  -  2) 

Choose  Ke  such  that  at  t  =  t0,  y  =  60  ~  w;  therefore, 


(6.23) 


K6  =  -  J60  (^«*  -  2)  . 


(6.24) 


To  obtain  the  initial  condition  for  t,  use  (6.10)  and  (6.18)  so  that 


1  =  t'o  —  2/3  Jecs  cos  y  —  Je2cs 


2  (15s2  —  14)(cos(2y  -  26)  -  cos(2w)) 


+  K2.  (6.25) 


24(5s2  -  4) 

It  is  evident  that  the  Je 2  terms  are  eliminated  when  y  =  6  -  u  and  6  =  60.  The 
result 


K2  =  2/3  Jecs  cos (60  —  u)  (6.26) 

gives  the  desired  t  =  t'o  at  t  =  to- 

For  H,  (6.10)  and  (6.19)  are  required.  In  addition  all  secular  terms  from 
fl2,  equation  (B.3),  App.  B,  are  needed  since  all  these  terms  are  of  order  J. 
Therefore, 


n 


fir,  —  J6c 


ce2J(75sA  —  102s2  +  56)(sin(2y  -  26)  -f  sin(2^)) 


24 (5s2  -  4)2 

4ceJ  sin  y  ce2  J2  cos(2u;)0  5ce2J2s2  cos(2u.')0 


15s2  -  12 


8 


ce7  J2  co$[2jj)6  ,  2  5 ce2J2s26  F>cJ2s26 

v  —  -  bcj2K2s26  - 


12 

ce2J20  c.J26 

--  -- 


24 


(6.2' 
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and 


4 

Kt  =  Jc60  +  -  Jet  sin(0o  —  w) 

3 


(6.28) 


8.  6  As  a  Function  of  Time.  (Initial  Conditions  for  r,  %j) 

To  complete  the  evaluation  of  the  remaining  constants  in  u,  it  is  nec¬ 
essary  to  give  0  as  a  function  of  time.  The  expression  for  the  perturbed  ^  then 
can  be  related  to  the  well-known  two-body  formula,  and  from  this  relationship,  an 
expression  for  r  and  ^  may  be  derived.  The  task  will  be  to  choose  Kg  and  Kg  so 
the  conditions  as  given  by  (5.27)  are  satisfied,  i.e., 

(t  ]  _  Q(!  ~  e2)  dr_(t  x  _  a(l  -  c2)e  sin(fl0  -  u) 

'  1  +  e  cos(^o  —  u)  d6  0  (1  +  e  cos(0o  —  u>))2 

Proceeding  in  this  manner,  the  formula  for  jfc  can  be  obtained  from  the 

relation 

dt  _  d<t>  dt  d<f>  r2  cos 2  6 
dd  dd  dtp  d$  h  cos  to 

Direct  substitution  of  the  derived  expression  for  equation  (5.21),  results  in 


At  t  =  tr 


r2  (  e2  cos(2y  -  26)  —  t2  cos(2u;)  e2s2  cos(2y  —  26) 

h  \  15  (5s2""—  4)  +  ~~  8  ~ 

c2  cos(2y  -  20)  es2  cos(y  +  26)  es2  cos(y  -  26) 

60  6  2 
4es2cos  y  4e  cos  y  s2cos(20)  2es2cos('^’  —  B~) 

3  3  2  3 

e2s2  cos(2u.')  e2cos(2u.')  s 2  \ 

8  +  60  +  ~2  ~  1  j 


(6.29) 
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The  above  condition  is  an  integral  from  the  two-body  problem.  If  (6.29)  is  evalu¬ 
ated  at  $  =  $0,  then  ||  may  be  expressed  as 


r<1+ £ 


where 


es2  cos(w  +  60)  2 es2  cos(w  -  0O)  4e  cos(u  -  60) 

K'° - 2 - + - 3 - J - 

es2cos(w  —  30o)  sJcos(250)  s2 

6  2  h  ~2  ~  1 

From  this  expression,  a  formula  for  h  results: 


so  that  from  (6.29) 


h  —  /i(l  +  J  K  io) 


Using 


P  =  53/  =  5S(1  +  2JKlo)  =  0(1  ~  '!)(1  +  2JK'° > 


a  new  formula  for  r,  which  includes  the  first  order  perturbation  effects,  may  be 


written  as 


r  =  P  =  _ o(l  -  e2) _ 

u  1  +  e  cos  y  -r  J (-2 Ki0  +  iti) 
where  Ui  is  (6.17). 


(6.30) 


The  formula  for  jjjj  is 

dr 


a(l  -  e2) e  sin  y 
(1  +  e  cos  y)2 

Ja{e“  -  1)  e2s2  sin(2y  4-  26) 
(l  +  e  cos  y)2  6 


5es2  sin(y  4-  26) 


-  eK7  sin(y  -  26) 
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„  .  ,  .  ,  lle2s2  sin(2y)  5e2sin(2y)  5es2sin(y) 

-  eKjsm(y  +  2u>) - 6  -  + - - 

.  .  e2$2sin(20)  s2sin(20) 

-  Kg  sin  y  +  2c  sin  y  +  Ks  cos  y  -f - — - — -  H - - - 

Z  o 


(6.31) 


Next,  evaluate  (6.30)  and  (6.31)  at  0  =  0O.  Keep  only  the  terms  to  order  J.  The 
result  is  two  simultaneous  equations  that  can  be  solved  for  Kg  and  K9: 


c2s2  cos(3a>  -  0o)  lle2s2  cos(3u;  -  30o) 

16  +  24 

5e2cos(3w  —  30o)  e2.s2  cos(3w  —  50o)  31es2  cos(2u>  —  20o) 

12  16  +  12 
7c  cos(2w  —  2 60)  3css  cos(2w  —  40q)  lle2s2  cos(w  +  0O) 


s2  cos(u;  +  0O)  llc2s2  cos(u>  —  0o)  7s2  cos(w  -  0O) 

- 4 - + - 5 - + - 2 - 


31e2  cos(w  —  0O) 


-  3cos(w  -  0O) 


17e2s2  cos(w  —  30o) 


7s2  cos(*' -  30o)  cs2  cos(2x')  5es2cos(20o)  13es2  7 e 


12  3 


(6.32) 


e2s2 sin(3~'  -  0C)  lle2s2 sin(3~’  -  30c)  5e2(3x’ —  30q) 

—  —  ^  — — 

e252  sin(3-c’ —  50o)  31e52  sin(2w  —  20o)  7e  sin(2^' -  20q) 

16  12  '  3 

3e,s2  sin(2^' -  40e)  7e2s2  sin(u;  +  0o)  s2  sin(w  +  0O) 

8  16  +  4 

67eJs2  sin(u;  -  0C)  7s2  sin(u>  -  0O)  29e2  sin(w  -  0O) 

24  2  12 

.  „  ,  1  Ie252  sinfu.’  -  30o)  7s2  sin(u?  -  30r) 

-  3  s  n  *  -  0r,  - - - - -  -4-  - i - - 

v  48  12 

C52sin(2^.')  3e52sin(20o) 

2  *  4  ' 


(6.33) 


Substitution  of  the  values  of  K%  and  K$  in  the  equation  for  u  insures 
that  the  remaining  initial  conditions  for  r  and  ^  are  satisfied. 

C.  THE  RESULTS 

1.  Preliminaries 

It  was  noted  that  throughout  the  calculations  in  this  chapter  the  term 
(s»‘n  tgr2)  appeare(j  jn  the  denominator  of  terms  in  Ui,  t2,  H2,  and  u2.  It  would 
appear  that  the  solution  is  not  valid  when  t0  =  ±sin-1  \Ja/ 5.  And  if  the  solution 
is  not  valid  near  these  “critical  inclinations” ,  then  one  should  question  the  validity 
of  the  underlying  perturbation  process. 

The  purpose  of  this  last  section  is  to  display  the  final  solutions  for  y,  t, 
fl,  and  u,  and  to  show  that  each  of  the  apparent  singularities  can  be  eliminated 
in  the  limit  as  as  t0  approaches  ±sin-1y/4/5.  It  is  a  remarkable  aspect  of  this 
analysis  that  every  term  that  appears  to  cause  a  solution  to  blow  up  is  exactly 
canceled  by  a  corresponding  term  that  is  “hidden”  in  the  solution. 

2.  First  Order  Solution  for  y 

By  use  of  trigonometric  identities,  (6.23)  may  be  rewritten  as 


y  =  6  —  w  -j-  J6  (  -s2  —  2 


-  Je 


2  (75s6  -  260s*  +  296s2  -  112)  /_  „  (b 


96  (| s2  -  2)' 


cos  (  2w  -  J6  (  ^s2  -  2^jj  sin  ( J6  -  2  j  j 


eV)te°5(2f  +-Jejy-  2)  +  J'H,  +  + 

30  (|s2  —  2)  V2  / 


(6.34) 


Equation  (6.34)  is  the  complete  first  order  solution  for  y.  If  the  initial 
inclination  is  such  that  ||s2  —  2  <  1CT3,  then  (6.34)  should  be  replaced  by  its 
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limit: 


e2J3  sin(2w)02  8eJ2  cos(u  -  0O)6  17e2J2  cos(2u)0 

y  _  -w  +  _  +  —  + 


+ 


e2J20  2 J20 

60 


+  2~Y  +  o  (j2,  J50,  J4ea  s 2  -  2)  + . . .) 


3.  First  Order  Solution  for  * 

As  was  done  with  the  expression  for  y,  i  is  rewritten  as: 

t  =  t0  +  J  c  s  ^e[cos(0o  —  w)  —  cos  y] 


e2(l5s2  —  14)  sin 

2 u  -  J6  | 

(!**-*) 

j  sin 

JO  | 

(I**-*)]) 

24  | 

f  s , 

52  -  2) 

1 

J 

+  U[J  ‘  +  Ja0  +  J*62  +  ...).  (6.35) 


Again,  as  in  y,  should  t0  be  near  the  critical  inclinations,  ||s2-2!  <  10  3, 
then  (6.35)  should  be  replaced  by: 

4  Je.  J20e2  sin  2u 

,  -  t0  +  — —  [cos(0o  -  u)  -  cos(0  -  w)j  4 - — - 

15  30 

+  O  (  J2,  J30,  J302  f- sin2  to  - 


4.  First  Order  Solution  for  U 
Rewriting  (6.27)  for  fl  results  in 


fl  =  fl0  —  J  6c 

_2  7(75s<  -  120s2  +  56) 

-f-  J  7) 

48  (fs2  -  2) 

cos  ^2x'  -  J6  ^s2  —  2^  sin  ^  J6  ^s2  ~  2 J  J 

4ceJ  sin  y  ce2J2  cos(2u,>)0  5ce2J2s2  cos(2u;)0 
3  15s2  -  12  8 
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ce2J2  cos(2u>)0 


+  5cJ2K2s28  - 


5  ce2J2s26  5  cJ2s26 


ce2J28  cJ26  ,  ,  , 

r - 1 - - - \-K4  +  0(J  ,J  6,J  6  +  . . .) 


(6.36) 


where 


Ki  =  -Jecs  cos(^o  —  w) 

O 

4 

K4  =  Jc80  +  -  Jce  sin(0o  —  w) . 

3 

If  j | s2  —  2j  <  10_s,  then 

4ceJ  sin(0  -  w)  cc*Js  sin(2a;)0i  bce2J2  cos(2w)0 

ll  =  Itn  —  Juc  —  -  + - 1 - 

3  6  12 


+  4cJ2K28- 


ce2J26  11  cJ26 


+  K4 


+  o  (V  +  J*e  +  j'e3  -  2)  + . . .) 


where 


K2  =  -7  Je  cos(0o  -  u) 

15 

Ki  =  \/l/5  J6 0  4 — i/l/5  Je  sin(0o  -  w)  . 
3 


5.  First  Order  Solution  for  u 


1  +  e  cos  y  - 


The  complete  first  order  solution  for  u  is: 
e2Js2  cos(2y  -  2$)  ,  0  ,  T  „  , 


+  2eJ  J KiS2  cos(2y  —  26) 


e2Js2  cos(2y  -  26)  5e  Js2  cos(y  4-  26)  lle2Js2cos(2y) 

+  - g - ^ -  +  - Y2 - eJKl  C°S'y  ~  20^ 


eJ K-j  cos(y  2 uj)  - 


5e2J  cos(2y)  e2Js2  cos(26)  Js2  cos(26 ) 


-  J/f&sinU'  -  0O)  -  JA'gcos(u;  -  0O)  -  2e2JKiS2  cos(2~’)  -  2 JKts2 


5  J?2  7e2J 


J  *  0(J2) 


7G 


where 


15s2  -  14 


48  (§s2  -  2) 

=  2 ecs  cos(60  —  u) 

=  8-u  +  je(^s'-  2) 


2  (75s6  -  260s*  +  296s2  -  112) 


(2w-je  (j*!-2))sm  (jb  (|**-2)) 


48  (|s2  —  2)  V  \2  JJ  \ 

eiJ28  cos(2w)  +  JQo  ^s2  -  2^  +  J20y 2  (see  equation  6.22  for  y2) 

—  15e2s*  —  (I4e2  —  2)s2 
48  (§s2  -  2) 

Equation  (6.32) 


- Equation  (6.33) 

Applying  the  condition  [|s2  —  2j  <  10~3,  the  expression  for  u  is 


u  =  1  — 


e2 J  cos(40  —  2u>)  e2J  cos(2u>)  eJ  cos(30  —  u) 


e2J  cos(2 9  -  2'J) 


JKg  s\n(60  —  w)  +  JK%  cos(0o  -  w) 


/e2J3  sin(2w)02  8eJ2  cos(u  —  6o)0 

-r  e  cos  - - — - - h - - - — 

v  15  15 

17e2J2  cos(2u;)0  e2J26  2  J2B  .  \ 

60  60  5  J 

e(8e2  —  8)  J26  sin(0  +  w)  t2J  cos(20)  2 J  cos(20) 

120  5  15 

_  gf _  MU 1  +  _  j  +  o  (j2,jse,j3e2  f-s2  -  2 

15  5  30  V 

Where  /f?  and  Kg  are  now: 

e‘  cos(3x-  -  60)  t2  cos(3^  -  30o)  e2  cos(3w  -  580) 


4t  cos(2u>  -  20o)  3e  cos(2u>  —  40o)  lie* cos(w  +  0O) 
15  10  20 

cos(u  +  $0)  89e*  cos(cj  —  0O)  cos(w  —  80) 

5  60  5 

17e*  cos(w  -  360)  7  cos(w  -  30o)  2e  cos(2u;) 


60 


15 


-  e  cos(20o) 


Ks 


e*sin(3w  -  0O)  ,  e*  sin(3u;  -  380)  i  e*sin(3w  —  50o) 
20  +  20  +  20 

4e  sin(2u  —  260)  3e  sin(2u;  —  40o)  7e*  sin(w  +  80) 

15  +  10  20 


sin(w  0O)  lie*  sin(cj  -  &o)  sin (u>  —  6o)  lie*  sin(w  -  380) 

+  5  +  60  +  5  +  60 

7  sin(u  -  30o)  2e  sin(2oj)  3e  sin(20o)  t2. 

+  - r~ - + - i - i - +  °yJ  ) 

lo  5  5 

The  results  obtained  thus  far  may  be  used  to  predict  the  orbit  of  a  satel¬ 
lite  for  up  to  1000  revolutions  when  a  valid  set  of  initial  conditions  are  provided. 
The  initial  conditions  will  usually  be  given  as  the  initial  displacement  vector  r 
and  the  initial  velocity  v,  whereas  the  solution  obtained  is  in  terms  of  the  orbital 
elements.  It  is  important  to  show  that  the  coordinates  and  velocity  components 
can  be  easily  recovered  from  the  orbital  elements.  Indeed,  Reference  4  contains 
the  criticism  that  Brenner  and  Latta  did  not  show  how  r  and  v  would  be  ob¬ 
tained.  and  that  ‘‘the  derived  elements  are  not  such  that  the  velocity  components 
are  readily  obtainable.’’  (Arsenault  et  al.,  Ref.  4:  Vol.  2,  p.  5). 

Again  referring  back  to  Figure  5.4  and  equations  (5.1)  -  (5.3).  the  po¬ 
sition  of  the  satellite  in  the  coordinate  system  may  be  derived  from  its  direction 


cosine;: 


r 


r  cos  f  cos  o  i  -  r  cos  t  sin  O  j  *  r  sin  6  k 


where  r  is  equation  (6.30)  and  i,j,k,  are  unit  vectors. 
Measuring  from  the  line  of  nodes  results  in 

r  =  r  cos  6  cos(d>  -  fl)i'  +  r  cos  6  sin(d>  —  ft)j'  +  r  sin  6  k 
where  i'  =  cos  fl  i  +  sin  fl  j 

j'  =  -  sin  H  i  +  cos  fl  j. 

Using  the  angle  relationships  from  equation  (5.1)  -  (5.3),  the  above  expression  is 
rewritten  as 


r  =  r(cos6  cos  fl  —  cos  t  sin  6  sin  fl) i 

+  r(c os  6  sin  fl  +  cos  i  sin  6  cos  fl)j  +  r  sin  i  sin  6  k  .  (6.38) 

The  velocity  is  found  by  differentiating  (6.38)  with  respect  to  6  and  noting 

dr  dr  dO 
dt  d6  dt 


The  result  is 

(  dr i .  dr 2 .  dr3\dd 

v  “  (d#1*  d0J ~  lek)  dt 

where 


dr  i 
d6 


.  .  dr  .  „  dfl  .  n 

-  cos  z  sin  fl—  sin  6  -  cos  t  cos  fl-— r  sin  6 
d6  d6 

.  di 

sin  t—  sin  Cl  r  sin  6  —  cos  fl  r  sin  6 
dO 

dr  n  .  dfl 
cos  fl  —  cos  6  —  sin  fl-— r  cos  6 
d6  d9 


—  cos  z  sin  fl  r  cos  6 


—  r 


.  .  dn  .  ^ 

-  cos  z  sin  fl  — - ■  sin  0 
dO 


sin  fl  sin  6 


dr  2 
~dS 
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di  dCl 

-  sin  t—  cos  n  sin  0  +  cos  D—  cos  0 
dO  dO 


+  cos  tcos 


dr 


n  cos  o^j 


+  —  (cos  * cos  n  sin  0  +  sin  fi  cos  6) 


drs  .  .  .  ,  .  .  .  .dr  .  .  di 

— —  =  r  sin  i  cos  0  +  sin  *  sin  0—  +  r  cos  t  sin  6 — . 
dO  dO  dO 


The  requirement  now  is  to  find  expressions  for  and  the  expression  for 

^  is  equation  (6.31). 


Differentiating  the  expression  for  i,  equation  (6.35),  and  the  equation 


for  fl,  equation  (6.35),  the  following  equations  for  ^  and  ^  are  obtained: 


di 

To 

2  Jecs  sin  y 

dQ 

4  T 

To 

ji 

—Jc  —  -  Jet  sin  y  . 
o 

(6.39) 

The  expression  for  ^  is  equation  (6.29).  By  expanding  this  expression, 
the  following  equation  for  ~  is  obtained: 


dO  h  (  e2  cos(2y  —  20)  -  e2  cos(2w)  e2s2  cos(2y  -  26) 

dt  ~  r2  1  _  \  15(5s2  -  4)  +  8 

e2cos(2y-20)  es2  cos(y  +  20)  es2  cos(y  —  26) 

60  6  2 
4es2cos  y  4e  cos  y  s2cos(20)  2es2cos(w  —  0O) 

3  3  2 

e2s2cos(2^)  e2cos(2w)  s2 
8  +  60  +  ~2 

For  j|s2  -  2  <  10~3 
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D.  VERIFICATION  OF  THE  RESULTS 
1.  Comparison  with  Brenner  and  Latta 

The  comparison  of  solutions  obtained  by  Brenner  and  Latta  must  be 
restricted  to  the  J,  Je,  Je 2,  and  J2  terms.  Brenner  and  Latta  limited  their 
analysis  in  J  to  these  terms  and  thus  avoided  the  secular  terms  that  arise  in  J2e2. 
No  mention  is  made  of  the  critical  inclination  in  their  work  because  ||s2  —  2|  does 
not  appear  as  a  divisor  in  the  terms  they  included.  Similarly,  this  present  analysis 
neglected  all  harmonics  in  the  potential  except  J2  while  Brenner  included  some  of 
the  terms  associated  with  the  JA  (or  Dx  by  their  notation)  harmonic. 

Before  comparing  solutions,  it  should  be  noted  that  there  is  some  dif¬ 
ference  between  the  two  works  in  the  notation  used,  this  analysis  preferring  the 
more  standard  Astrodynamic  notation.  Brenner  used  the  co-latitude  angle  (0)  as 
a  spherical  coordinate  while  the  latitude  (<5)  was  used  here.  In  addition,  Bren¬ 
ner  used  M  as  the  independent  variable  where  M  +  n/2  defines  the  angle  from 
the  ascending  node  to  the  satellite.  The  independent  variable  6 ,  measured  from 
the  ascending  node  to  the  satellite,  was  used  here.  Finally,  Brenner  defined  the 
rotation  of  the  line  of  nodes  in  an  opposite  sense  to  that  done  here.  In  summary, 

6  — >  7r/2  -  6  (where  6  is  co  —  latitude) 

0  — ♦  M  +  7r/2  (where  6  is  the  polar  angle) 

n  —  -n 
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There  w;U  also  be  a  difference  in  the  constant  terms.  Brenner  chose 
initial  conditions  that  would  make  the  analytical  solutions  as  simple  as  possible. 
This  analysis  adopted  the  more  general  initial  conditions  that  were  outlined  in 
Chapter  5. 

Now  restricting  the  comparison  to  the  terms  that  are  left,  Brenner’s 
solution  for  i  and  fi  are  respectively: 


i  —  — Jecs  cos  y  +  . . 
3 


n  =  JMc- 


J2Mc 


- J2Mcs 2  -I-  -Jtc  sin  y  + 

3  3 


5  J2 


cs 


12 


sin(2m)  + 


A  check  of  equations  (6.35)  and  (6.36)  shows  there  is  an  exact  match  of  terms  of 
order  J.  The  final  term  in  Brenner’s  solution  for  fl  is  a  J2  term  and  can  be  found 
in  the  solution  for  fl2,  equation  (B.4).  Next,  Brenner’s  abbreviated  solution  for 
is 

5es2  s2  5 

- cos(y  +  2M)  -\ - cos(2Af)  -  -s2  -f  1  4-  . . . 

24  6  2 

These  same  terms  are  duplicated  in  equation  (6.17).  Finally,  Brenner’s  solution 

for  yi  is  |s2  —  2  which  is  exactly  the  expression  obtained  in  this  analysis. 

2.  Comparison  with  an  Independent  Analysis  of  Polar  and 
Equatorial  Orbits 

Appendix  A  contains  an  independent  analysis  of  the  polar  and  equatorial 
orbits.  Since  for  the  polar  case  there  is  no  variation  of  inclination,  no  rotation  of 
the  line  of  nodes,  nor  dependence  of  the  motion  on  the  longitude  <j> ,  new  equations 
must  be  derived.  The  equations  cannot  be  solved  in  terms  of  the  angles  i  and  H. 
Instead,  the  equations  are  solved  in  terms  of  the  variable  A2  wrhich  is  related  to 

di 

d6  ' 

For  the  equatorial  case  the  inclination  remains  constant.  The  line  of 
nodes  is  undefined  since  the  satellite  remains  in  the  equatorial  plane;  however,  the 
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angle  0  does  change  with  time.  Instead  of  using  the  polar  angle  6  to  measure  the 
angular  displacement  of  the  satellite,  the  angle  <f>,  measured  from  the  fixed  axis  7, 
is  used.  A  new  dependent  coordinate  is  defined  as 

Y  =  4>  -  K  +  JYrf  +  J2Y,<i> . 

As  Appendix  A  demonstrates,  there  is  exact  agreement  between  the 
special  polar  case  and  the  general  case  for  i'o  =  90°.  In  addition,  when  the  appro¬ 
priate  change  of  variables  is  made  in  the  equatorial  case,  it  agrees  exactly  with 
the  general  case  at  t0  =  0°. 
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VII.  CONCLUSIONS 
AND  RECOMMENDATIONS 

The  results  achieved  in  Chapter  VI  represent  a  unique  solution  to  the  problem 
of  a  satellite  in  orbit  about  an  oblate  planet.  In  fact,  prior  to  this  present  work 
the  problem  had  not  been  solved  in  this  representation  [Ref.  18,  p.  340].  In 
contrast  to  the  universal  solution  presented  here,  all  other  methods  require  a 
reformulation  of  the  problem  in  alternate  variables  in  order  to  achieve  solutions 
at  various  singularities,  the  most  prevalent  being  the  critical  inclination. 

The  results  support  the  theory  that  at  the  critical  inclination  there  are  no 
discernable  features  in  the  perturbations  of  the  elements  that  distinguish  it  from 
any  other.  This  conclusion  was  reached  by  Lubowe  [Ref.  38],  is  shared  by  Taff 
[Ref.  18]  and  has  been  corroborated  by  the  physical  data.  However,  all  analyti¬ 
cal  solutions  arrived  at  via  canonic  transformations  predict  perturbations  in  the 
vicinity  of  the  critical  to  be  25  times  greater  than  those  away  from  it. 

The  perturbation  method  used  in  the  analysis  embodied  the  principles  out¬ 
lined  in  Chapter  I.  First,  the  resulting  solutions  are  significantly  mere  accurate 
than  the  two-body  solution.  For  the  appropriate  orbits,  the  relative  error  of  two- 
body  solution  is  1000  times  greater  than  the  solution  obtained  here.  Second  the 
solution  was  obtained  in  parameters  closely  related  to  the  classical  orbital  elements 
and  in  Cartesian  coordinates;  no  transformation  to  an  alternate  non-physical  set 
of  elements  w-as  required.  Therefore  the  physical  effects  are  easily  distinguishable 
throughout  the  analysis.  Finally,  as  has  been  noted  the  solution  is  valid  for  all 
orbital  parameters. 
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While  the  perturbation  method  was  similar  to  the  Method  of  Strained  Co¬ 
ordinates  several  novel  techniques  for  dealing  with  secular  terms  and  apparent 
singularities  were  introduced.  The  verification  of  these  techniques  was  possible 
due  to  the  extensive  use  of  the  symbolic  computer  program  MACSYMA  which 
allowed  for  the  investigation  of  higher  order  formulas. 

The  analysis  suggests  several  areas  for  further  research.  The  areas  include: 

•  The  addition  of  the  higher  harmonics  of  the  gravitational  potential,  e.g.,  Js, 
J4,  etc.,  to  the  problem. 

•  The  investigation  into  the  feasibility  of  including  non-conservative  perturb¬ 
ing  forces  such  as  drag  or  solar  radiation  pressure. 

•  Numerically  integrating  the  equations  of  motion  subject  to  the  initial  condi¬ 
tions  used  here  in  order  to  check  the  analytic  results. 

•  An  in-depth  analysis  of  the  use  of  canonic  transformations  in  satellite  orbit 
prediction  to  include  a  comparison  with  the  method  used  here,  the  purpose 
being  to  determine  if  the  critical  inclination  singularity  is  an  artifact  of  the 
transformation  process. 
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APPENDIX  A 


THE  POLAR  AND  EQUATORIAL  ORBITS 


A.  THE  POLAR  ORBIT 

Reference  is  made  to  the  general  equations  of  motion,  equations  (5.13)  - 
(5.15).  Now  in  the  case  of  a  polar  orbit,  the  longitude  <f>  is  a  constant.  Therefore, 
the  equations  of  motion  are  reduced  to 


dl 

dt 2 


r  (ddV  -GM  „  T  R2  /I  3  .  ,A  , 

?-rUJ  *)  <A1) 


,d0 
dt  (  dt 


GM  R2  .  .A. 
—  Ji - ; —  sin(20) 


Let  J  =  |  (^2^)  and  u  —  p/r 


where  p  =  h2 /GM. 


(A.2) 


Define  a  function  A  such  that 


so  that 


2M  t  a 
r  —  -  h  A 
dt 

d  hu2 A  d 


dt  p2  dO 

Applying  (A. 3)  to  (A.2)  the  independent  variable  can  be  changed  from 
Equation  (A.2)  may  be  rewritten  as 

dA2 


de 


Similarly  equation  (A.l)  becomes 


-  — 2Ju  sin(20). 


(A- 3) 
t  to  e. 

(A.4) 


d2u  1-Ju[^sin(2  + 

d6  A2 


(A.5) 
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To  solve  (A.4)  and  (A. 5)  a  solution  of  the  form 


A2  =  1  +  J Ei  ■+■  J7 Ei  +  . . . 

(A.6) 

u  =  1  +  e  cos  y  +  Jui  +  J2U2  +  .  • . 

(A.7) 

where  y  =  6  —  u  +  J 6yx  +  J76y2  +  . . . 

(A.8) 

is  assumed. 

Substitution  of  equations  (A.6)  and  (A. 7)  into  (A.4)  and  equating  coefficients  of 
order  J  results  in: 

dE 

— — -  =  —e  sin(y  +  26)  +  e  sin(y  -  26)  -  2sin(20).  (A.9) 

du 

A  solution  of  (A.9),  ignoring  terms  of  order  J2  or  higher,  is 

„  e  cos (y  +  26)  .  .  .. 

Ex  =  - - -  +  e  cos(y  —  26)  +  cos(20) 

+  Kxe7  cos(2y  —  26).  (A. 10) 

Substituting  (A. 10),  (A.6),  and  (A. 7)  into  (A. 5)  yields  to  order  J 

dux  5e2(2y  +  2 6)  ,  . 

—  +  «i  =  2ecosyyi  +  - - - e  Kx  cos(2y  -  26) 

aU  o 

e2cos(2y-20)  5e  cos(y  +  26)  e2cos(2y)  . 

8  3  4 

3e2cos(2 6)  cos(20)  e2  1 

4  ^  2  7  ~  2' 

The  cos  y  term  would  produce  secular  terms  in  the  s  .ion  to  Uj,  therefore  yx  is 
chosen  as  ^  to  make  the  coefficient  of  cos  y  zero. 
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The  resulting  equation  is 


dui  ,  5e2  cos(2y  4-  20)  Xvr  e2cos(2y-20) 

_-  +  Ul  =  - ^ - -  -  t'Kx  cos(2y  -  20)  + - K-~ - >- 


5e  cos(y  +  20)  e2cos(2y)  3e*cos(20)  cos(20) 

+  3  ~  4  +  4  "*  2 


e2  1 

7  ~  2‘ 


It  is  assumed  that  the  solution  to  (A. 11)  has  the  following  form: 


ui  =  a0  +  oi  cos(2y  -f-  20)  +  a2  cos(2y  —  26) 


(A.ll) 


+  a3cos(y  +  26)  +  a5  cos(20). 


(A.12) 


Equation  (A.12)  is  substituted  into  (A.ll).  The  unknown  coefficients  may 
then  be  solved  for  by  equating  the  coefficients  which  have  the  same  harmonics. 
The  results  are: 


ao  =  --  -  - 


«1 


g2  =  - 


a3  =  - 


e*  1 
4  2 

24 

(8/fi  -  l)e2 


8 


5e 

24 


a<  =  — 


«5 


e 

12 

e^  l 

4  6 


Substitution  of  these  coefficients  into  (A.12)  results  in 

e2  cos(2y  26)  -  e2  cos(2y  -  20) 

= - — —  -  e2KiCos(2y  -  29)  +  —  1 


24 


8 
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5e  cos(y  +  20)  e2cos(2y)  e*cos(20)  cos(20) 

~  24  +  12  4  6 

f!_l 

4  2  ' 

The  following  complementary  solutions  must  be  added  to  the  particular  so¬ 
lution  for  u  to  ensure  that  any  secular  terms  in  the  second  order  solution  can  be 
eliminated 

—Ks  cos(3y  —  20)  —  Kj  cos(y  -  26) . 


The  first  order  solution  for  u  now  has  the  following  form: 

,  t  Ty  ir.  „a\  cos(2y  +  26) 

u  -  1  +  e  cos  y  —  J  Ks  cos(3y  —  26) - ' 


24 


e2  J (8 Ki  —  1)  cos(2y  —  26)  5 eJ  cos(y  +  26) 


8 


24 


—  J  Kj  cos(y  —  26) 


e2J  cos(2y)  „ 

+  - — - -  +  J  Kg  sin  y  +  J  K»  cos  y 

(3e2  -f  2)  J  cos(20)  e2J  J 
12  ~4  ~  2 ' 


(A.  13) 


The  expressions  for  Ex  and  are  substituted  into  (A.6)  and  (A  7)  respec¬ 
tively.  (A. 6)  and  (A. 7)  are  then  substituted  into  (A. 4).  The  result  is 
^  -  -Ks  sm(3„-^)  +  e'Sin(2!'  +  'll,)  c3  sin(2y  +  20) 


do 


24 


12 


, e2sin(2y  —  26)  , 

+  e2hi  sin(2y  -  26)  -} - — - e2Kx  sin(2y  -  46) 

A  4 

e2 sin(2y  -  46)  5esin(y  +  40)  ,  e  sin(y  +  26) 

+  - g - +  - ^ - Ki  sin(y  +  26)  +  - - - 

+  K9  cos(y  +  26)  +  Agsin{y  -  26)  f  C  S'n^ — —  -  Jf9cos(y  -  26) 


K7  sin(y  -  46)  +  Ks  sin(3y)  +  t2Kx  sin(2y)  - 

„  .  5e  sin  y  e2sin(40Gx)  sin(40) 
A7siny - — —  +  - - - -  + 


e2  sin(2y) 


24 


eJsin(20) 


sin  (20) . 


89 


Integrating  this  expression  for  Ej  would  yield  a  secular  term  in  the  harmonic 
cos(2y  -  20).  Therefore  the  constant  Ki  is  chosen  as  —  ^  in  order  to  eliminate  the 
problem  term.  The  equation  is  now 


dE2 

d$ 


=  -K,  sin(3y  -  40)  +  +  46)  -  +  28> 

v  '  24  12 

5e2  sin(2y  -  40)  5e  sinfy  +  46) 

+  - ^ - -+ - g - ^-/f8sin(y  +  20) 

e  sin(y  4-  20)  .  e  sinfy  —  20) 

~h  - — - f-  i£g  cos^y  +  20)  +  /fg  sin(y  —  20)  4 - — - 

-  Kg  cos (y  -  20)  -  AT7  sin(y  -  40)  +  Kh  sin(3y)  - 


4-  K 7  sin  y 
e2sin(2  0) 


5e  sin  y  e2sin(40)  ^  sin(40) 


24 


4-  sin(20). 


Integrating  the  above  expression  yields 

£,  =  -  g,  cos(3y  -  4«)  -  C°S(2y  +  4>)  - 

144  48 

5e2  cos(2y  -  40)  e  cos(y  +  40)  Kg  sin(y  4-  20) 

48  24  +  3 

K$  cos(y  -J-  20)  e  cos(y  4-  20) 

_  — 


4  Kg  sin(y  -  20) 

K  cosf  -  20]  *  e  cos(^  ~  26>)  __  #7Cos(y  -  4x)  _  cos(3y) 

e2cos(2y)  5e  cos  y  e2cos(40)  cos(40) 

- K7  cos  y  4 - — - — - - — - 

C  tv  *  .  i 


24 


24 


e2cos(20)  cos(20) 

-  j 


(A. 14) 


Continuing  the  procedure  results  in  the  following  second  order  expression 

d2u2  e/C5  cos(4y  -  20)  e4cos(4y-20) 

— —  u2  =  2e  cos  yyo  - - - - - - - — - - 

d62  2  96 

eA'j  cos(4y  -  40)  5e4  cos(4y  -  40)  3e3  cos(3y  —  40) 


576 


16 


90 


3es  cos(3y  +  20) 


—  2 Ki  cos(3y  —  26)  — 


e*cos(3y  —  26) 


35e2  cos(2y  +  46)  5eK9  sin(2y  +  26)  5cK8  cos(2y  +  26) 
32  +  4  +  4 

5e2  cos(2y  +  26)  eK9  sin(2y  -  26) 

12  +  4 

eKgcos(2y  -  20)  tK7  cos(2y  -  26)  tKg  cos(2y  -  26) 

4  +  ~2  +  2 

e*  cos(2y  —  26)  7e2  cos(2y  —  26) 

48  +  24 

3eK7  cos(2y  —  46)  3eKg  cos(2y  —  46)  e4cos(2y  —  46) 
~  4  4  +  32 

e2cos(2y-40)  7escos(y  +  40)  7c  cos(y  +  40) 

8  8  4 

5Kg  sin(y  +  26)  5A"gcos(y  +  26) 

3  +  3 


e3cos(y  +  20)  13e  cos(y  +  2 6) 


+  2 K7  cos(y  -  20) 


16  36 

e3cos(y  — 20)  e  cos(y  -  20)  5K7  cos(y  —  40) 

16  6  3 

55e3  cos(y  -  40)  5ecos(y-4 0)  5eKs  cos(4y) 


144 

12 

4 

5e4  cos(4y) 

5 K*,  cos(3y) 

13e3  cos(3y) 

e/Cs  sin(2y) 

192 

3 

144 

2 

eK$  cos(2y) 

3 eK7  cos(2y) 

3eK$  cos(2y) 

2 

4 

4 

e4  cos(2 y) 

23e2  cos(2y) 

25e3  cos  y 

32 

32 

48 

17e  cos  y 

5eK/cos(46) 

4 

5e4cos(40)  65e2cos(40) 

24 

192 

32 

5cos(40) 

8 

3eKg  cos(20) 

2 

eK7  cos(2 0) 

2 

e4cos(20) 

37c2  cos(20) 

cos(20)  cXg 

96 

48 

3  2 

eK7  5e4  167e2  1 

4  576  28&  ^  6' 


(A. 15) 


9] 


It  is  noted  that  three  harmonics  in  (A. 15)  would  p.^duce  secular  terms  in 
the  solution  to  U2.  The  harmonics  are  cos(3y  —  20),  cos(y  —  20),  and  cos  y.  These 
harmonics  may  be  eliminated  by  choosing 


Ks 

K7 

V2 


96 

(3es  -  8e) 
96 

25e2  +  34 
96  ' 


The  first  order  solution  for  (A .4)  is 


A2  =  1  +  J  (cos(20)  +  e  cos(y  -  20) 

e  cos(y  +  20)  e2cos(2y  — 20) 


(A.16) 


u  may  now  be  rewritten  as: 


e2J  cos(2y  +  20)  5 e2J  cos(2y  —  20)  5eJ  cos(y  +  20) 

u  —  1  +  e  cos  y - - — — - -  +  - - - 

24  24  24 

e3J  cos(y  -  20)  tJ  cos (y  -  20)  e2J  cos(2y)  e2J  cos(20) 

+  25  12  +  12  5 


J  cos(20)  „  .  .  .  .  T  r  .  .  .  e2J  J 

- J  K9  sin(w  -  0O)  +  J  K$  cos(w  -  0O)  -  - - — 

6  4  2 


(A. 17) 


,  J0  Je3  sin(2y  —  20)  T?/,(25e2  + 

where  y  =  0-  w+  ~  +  - — - -  +  J20^ - — r 

2  48  96 


Equations  (A.16)  and  (A. 17)  satisfy  the  equations  of  motion  to  first  order, 
and  the  equations  produce  no  secular  terms  to  second  order. 

A  comparison  may  now  be  made  between  the  polar  solution  derived  in  this 
appendix  and  the  general  solutions  for  y  and  u,  equations  (6.22)  and  (6.37)  respec¬ 
tively.  To  make  the  comparison  easier,  the  relatively  complex  initial  conditions 
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(equation  5.27)  are  discarded.  Instead,  following  Brenner  and  Latta,  the  initial 
cond  Mons  are  selected  in  order  to  make  the  expressions  as  simple  as  possible. 
In  <.  iition,  constant  terms  like  sin(2cj)  and  cos(2u>)  are  dropped.  It  will  be  re¬ 
membered  that  these  constants  were  added  to  the  general  solution  to  prevent 
singularities  at  the  critical  inclination.  It  is  a  simple  exercise  to  add  the  cos(2w) 
term  to  the  polar  case,  specifically  to  equation  A. 10,  and  then  show  that  identical 
results  are  achieved  with  the  general  case  at  i0  =  90°.  In  fact,  this  has  been  done. 
However,  to  display  these  terms  adds  nothing  to  the  comparison  and  only  serves 
to  make  the  analysis  more  tedious. 

Checking  equations  (6.35)  and  (6.36),  it  is  seen  for  the  polar  case,  cos  i0  = 
c  =  0,  there  results 

i  =  i0  and  fi  =  H0 . 

This  agrees  with  orbital  theory  that  for  the  polar  case  there  is  no  variation  of  the 
inclination  nor  does  the  line  of  nodes  rotate. 

Setting  sin  i0  =  s  =  1  and  disregarding  the  appropriate  constants,  equations 
(6.37)  and  (6.22)  are 

j  ___  e2J  cos(2y  29)  be7J  cos(2y  -  26)  5 eJ  cos (y  —  26) 

e  cos  y  —  -  -  —  — 

e3J  cos(y  -  26)  eJ  cos(y  -  26)  e2J  cos(2y)  e7J  cos{26) 

24  12  ~  12  4 

J  cos(2^}  e1 J  J 

-  - - - J  K9  sin(u;  -  90)  +  J  Kg  cos(u;  -  60) - 

o  4  2 

and 

y  =  t-u+i£  +  -  m)  +  ±  30 

2  48  96 

These  results  agree  exactly  with  the  (A. 17)  and  (A. 18). 
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B.  THE  EQUATORIAL  ORBIT 

Reference  is  again  made  to  the  general  equations  of  motion,  (5.13)  -  (5.15). 
An  orbit  that  remains  precisely  in  the  equatorial  plane  is  possible  if  only  even 
harmonics  of  the  potential  are  considered.  Now  6  =  0°,  j~  =  0,  =  0,  and  the 

equations  reduce  to 

d2r  d<f>  dV  , 

^T-r^7  =  --3r  (A-18) 


d2r 

d<t> 

dt2 

~r  dt 

d 

(  j  dd 

dt 

As  in  the  general  case  and  the  polar  case  the  independent  variable  is  changed 
by  use  of  the  integral  obtained  from  (A.  19) 


2  d<i>  - 

r  —  =  constant  =  h 

dt 

d  u2  -  d 

dt  p 2  dp 


(A. 20) 


Substitution  of  (A. 20)  into  (A. 19),  results  in 

d2u  j 

—  =  1-  Ju2 

do~ 


(A.21) 


where  :  J  = 


IfjSL 

2  [Jl  p> 


To  solve  (A.21).  assume  a  solution  of  the  form 

u  =  1  —  e  cos  }’  +  Jui  +  J:u2 
where  Y  =  6  -  u‘  —  J <t>Yi  J2oY2. 


(A.22) 


Substituting  (A.22)  into  (A.21)  and  equating  coefficients  of  order  J  results  in 
»2 

^  i  r  i  r  ir  nr.  .  .  t  *  .2  2  r  /  k 


Ui  =  1  -  2 J  Yxe  cos  Y'  —  2 J e  cos  Y  —  e  cos  Y. 


(A.23) 


9-5 


The  presence  of  the  term  cos  Y  will  cause  secular  terms  in  the  solution  for  ui, 


therefore  choose  Yx  =  —1.  Equation  (A.23)  is  now 

<f2ui  e2  e2cos2, 

— +  tti  =  1  +  -  +  2 

Assume  a  particular  solution  to  the  above  equation 


(A.24) 


Ui  =  a0  +  Oi  cos  2 Y. 


(A. 25) 


Substitution  of  (A. 25)  into  (A.24)  results  in  the  following  solution  for  Ui 

«i  =  1  +  j  ~  ge2cos(2 y). 


(A. 26) 


The  approximation  for  u  is: 


e2  1 


u  =  1  +  e  cos  Y  +  J  (  1  +  — - e2  cos  (2Y) 

\  2  6 


(A. 27) 


The  solution  must  be  checked  to  ensure  it  causes  no  secular  terms  to  second  order. 
Substituting  (A. 26)  and  (A. 22)  into  (A. 21)  results  in 
d2u2  _  escos(2y)  ,  , 


+  Uj  —  —  ■ 


-f  e2  cos (2 Y)  +  2e  cos  Y  Y2 


5e3cos  Y  ,,  , 

+  - +  3e  cos  Y  +  e  +2. 

6 

Note  that  the  terms  with  cos  Y  will  cause  secular  terms  in  u2.  Therefore,  Y2  is 


chosen  as: 


5e^  _  3 

'72  ~~  2' 


The  first  order  solution  for  u  is  equation  (A. 27) 


where  Y  =  <M^2+2/' 


5e2  3 


(A. 28) 


Equation  (A. 27)  is  the  correct  solution  in  terms  of  the  variable  <t>.  To  compare 
the  solution  found  here  to  the  general  solution  for  the  case  (5  =  0°  will  require  that 


(A. 27)  and  (A. 28)  be  modified  such  that  they  are  in  terms  of  6. 
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Referring  back  to  Figure  (5.1),  it  is  seen  that  <f>  is  measured  from  the  fixed 


axis  T  such  that 

4>  =  n/2  +  6  +  n.  (A. 29) 

Now,  from  equation  (6.27),  the  first  order  solution  for  H  with  i0  =  0°  (sin  i0  = 

s  =  0  and  cos  i0  =  c  =  1)  is 

Ue1  4  J3e20  J36 

n  =  fi0  -  JO  H - —  sin(2y  -  20)  -  -Je  sin  y - - - b  — -.  (A.30) 

48  3  6  2 

(Note:  As  in  the  polar  case,  the  constant  terms  cos(2u>)  and  sin(2w)  have  been 

dropped  in  the  general  case  for  comparison  purposes.  As  explained  previously, 

eliminating  these  terms  does  not  affect  the  validity  of  the  comparison.) 

Substitution  of  (A.30)  and  (A. 29)  into  (A. 28)  results  in 


7Je3  7  4 

Y  —  6-u-2J0-\ - —  sin(2y  -  26) - -  J2e20  -  -  Je  sin  y. 

48  12  3 

Now,  from  (6.22),  the  solution  for  y  with  i0  =  0°  is 

7  7 

y  —  6  -  u  -  2J0 - J2e20  H — -Je2  sin(2y  -  2z). 

y  12  48  v  ' 

Using  equation  (A. 32),  equation  (A. 31)  may  be  written  as 

4 

Y  =  y - Je  sin  y. 

3 

And  therefore  by  use  of  the  Taylor  series  expansion 

4 

cos(F)  =  cos  y  +  -  Je  sin 2  y. 

O 

Equation  (A. 27)  may  now  be  written  in  terms  of  y(6)'- 


(A.31) 


(A. 32) 


u  =  1  +  e  cos  y  -r  J 


(  7  2  5 

1  +  7 

\  6  6 


e  cos  2  y  )  . 


(A. 33) 


Reference  now  is  made  to  equation  (6.37),  the  general  solution  for  u.  Ry 
letting  sin  »0  =  s  =  0  in  this  equation,  it  is  easily  seen  that  the  general  expression 
becomes  (A. 33). 
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APPENDIX  B 


2nd  ORDER  ANALYSIS 


Substitution  of  (6.2)  -  (6.5),  (6.10),  (6.16),  and  (6.17)  into  (5.23)  results  in 


.  dCl  .  (Pt2  .  .  ,d?Cl2  .  dt2 

—  2sm  i— -  sin  6 - -smH  sin  i  — —  cos  0  —  2—  cos  0 

de  d6 2  dO 2  dd 


=  cesKiss  sin(3y  —  6)  —  cK$s  sin(3y  —  0)  + 


ce5Ki  s  sin(3y  -  0) 


teWjmlSy-M)  +  cKis s.n(3y  _  3()  + 

3  6 

7 ce2ss  sin(2y  4-  30 )  143ceJss  sin(2y  +  0)  29 ce2s  sin(2y  4-  0) 

24  +  72  18 


4-  5ce2K3s3  sin(2y  —  0)  +  l0ce2Kxs 3  sin(2y  —  0)  - 


15ce2s3  sin(2y  —  0) 


-  4ce2/lT3Ssin(2y  -  0)  —  4ce2K\s  sin(2y  -  0)  4- 


7ee2s  sin(2y  —  0) 


Jis  to.,  o a\  ,  5ce2s3  sin(2y  30)  2 


5ce2K3ss  sin(2y  -  35)  4- 


4-  4ce2K3s  sin(2y  -  30) 


_  35ces3sin(y  +  35) 


-  4ce2K1ssin(2y  -  30)  - 


4-  ce  K7s  sin(y  4-  0  4-  2.j) 


2cesKiS 3  cos(2w)  sin(y  +  5)  —  ce3KiS  cos(2w)  sin(y  4-  5) 


3  •  /  ,  ,  5ces2  sin(y  4-  5) 


-  2ceK2s3  sin(y  4  0) 


4-  2cK^s  sin(y  4-  0) 


,,  .  .  ..  4ce  s  sin(y  4-  5)  „  , 

-  ceK2s  sin(y  4-  0) - - — — - 2cKgS  cos(y  4-  0) 

O 


-  ceK7s  sin(y  —  0  2u>)  + 


10ce3if1s3  cos(2w)  sin(y  —  0) 


cesKiS  cos(2w)  sin(y  -  0)  !0ceK2ss  sin(y  -  0)  5cesKiS3  sin(y  -  6) 

3  +  3  1  3 


41ces3sin(y  -  0) 


2c/f?ssin(y  -  0)  -  ceK7s  sin(y  -  0) 


ceK2s  sin(y  -  0)  ce3K\S  sin(y  -  0)  8cessin(y-5) 
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+  2 cK9s  cos(y  —  0)  —  eesXiS8  sin(y  -  30)  + 


9cesssin(y  -  30) 


-f  ceKjs  sin(y  —  30)  — 


cesK1ss\n(y  -  30)  5ce2ss sin(30)  5esssin(30) 


—  10ccjKi5S  cos(2u?)  sin  0  —  10 cK2ss  sin  0  + 
lOesssin0  ee2sin  0  .  . 

.  _j - - - cs  Sln  q 


5 ce2ss  sin  0 
12 


(B.l) 


3  3 

Choose  Ki  =  and  Ks  =  7S^2a°^t56  to  eliminate  the  secular  terms 

sin(2y  —  0)  and  sin(2y  —  30)  in  equation  (B.l),  so  now 


.  .df)2  .  „  d2x2  .  .  .d2Q2  dt2 

—  2  sin  x  — —  sin  0 - —  sin  0  +  sin  t  — ■  —  cos  0  —  2—  cos  0 

d0  d02  d02  d0 

39cess  sin(3y  —  0)  ee3s3  sin(3y  —  0)  .  , 

1800s2  -  1440  8  v  ' 

llcess  sin(3y  —  0)  35cess  sin(3y  —  30)  5ee3s3sin(3y  —  30) 


240 

+  cK5s  sin(3y  —  30)  + 


1800s2  -  1440  24 

7ce3s  sin(3y  -  30)  7ce2ss  sin(2y  +  30) 


+ 


144  24 

143ce2s3 sin(2y  +  0)  29ce2s  sin(2y  +  0)  35cess  sin(y  +  30) 


72 


18 


24 


+  ceK7s  sin(y  +  0  +  2w)  + 


78cess  cos(2u>)  sin(y  +  0) 


1800s2  -  1440 
ce3s3  cos(2w)  sin(y  +  0)  llce3s  cos(2w)  sin(y  +  0) 

4  120 

,  „>  5ces3  sinfy  +  0)  ,,  .  . 

—  2ceK2s3  sin(y  4-  0)  4 - - - -  +  2cKgS  sin(y  4-  0) 

8 

.  ,  ..  4ces  sin(y  +  0)  „  T.  , 

-  ceK2s  sm(y  4-  0) - 2cKgs  cos(y  4-  0) 


—  ceK7s  sin(y  -  0  4-  2w)  - 


3 

7 0ce3s  cos(2u/)  sin(y  —  0) 


4- 


1800s2  -  1440 
5cesss  cos(2w)  sin(y  -  0)  7 cess  cos(2u>)  sin(y  -  0) 

12  72 

35ce3s  sin(y  -  0)  lOcei^s3  sin(y  -  0)  5cess3 sin(y  -  0) 


1800s2  -  1440 


24 
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—  2cKgs  sin(y  —  8)  -  ceKjS  sin(y  —  0) 


41eesssin(y  -  0) 

12 

ceK^s  sin(y  -  8)  7 ce3s  sin(y  -  8)  +  8ees  sin(y  -  8)  +  2cK~s  r  '  8) 

3  144  3  »  [y 

39cess  sin(y  -  38)  cess3  sin(y  —  30)  9cesssin(y  -  38) 

+  18005*  -  1440  8  +  4 


+  eeKVs  sin(y  —  30)  — 


llee3s  sin(y  -  30)  5ce*s3  sin(30)  5csssin(30) 


240cels  cos(2u>)  sin  8  5ce*s3  cos(2w)  sin  8  ee*s  cos(2w)  sin  8 

+  1800s*  -  1440  4  +  6 

„  ,,  .  .  .  5ceiss  sin  8  10cs3  sin  8  ce*s  sin  8  .  .  . 

—  10cK2s3  sin  6  + - H - cs  sin  8 .  (B.2) 

12  3  3  v  ' 

Solve  for  x2  and  02  using  the  technique  described  in  Chapter  6,  Section  5.  The 

second  order  expressions  for  t2  and  fl2  are  respectively 

_  llce3s  cos(3y  —  28)  ce3s3  cos(3y  —  20)  23cess  cos(3y  —  20) 

*2  “  '  900s*  -  720  6  +  360 

2 cK5s  cos(3y  -  20)  25ee*s3  cos(2y  +  20)  29ce*s  cos(2y  +  20) 

3  96  +  144 

llce3s  cos(y  -  20)  2 cKss  sin  y  22cess  cos(2u>)  cos  y 

+  900s*  -  720  ~3  900s*  -  720 

ce3s3  cos(2u>)  cos  y  23cess  cos(2u>)  cos  y  8 ceK2s3  cos  y 

+  j  180  +  3 

ce3s3  cos  y  85ces3cos  y  2 cK&s  cos  y  2 cK7s  cos  y 

_  _  __  _  3 

2 ceKiS  cos  y  23ce3s  cos  y  20ces  cos  y  . 

3  360  9  v  ' 


n2  = 


75 s*  -  60 

4cKj  sin(3y  -  20 

4  30 

')  17ce*s*  sin(2y  +  20)  29ce*sin(2y 

3 

7 ces*  sin(y  +  20) 

72  '  144 

2ce3  sin(y  -  20)  ce3s*  sin(y  -  20) 

36 

3ces*  sin(y  —  20) 

75s*  —  60  12 

2c/f7  sin(y  -  20)  lice3 sin(y  -  20) 
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4ccscos(2cj)  sin  y  ce3s2  cos (2u>)  sin  y  eescos(2w)  sin  y 
75s2  -  60  +  2  15 


4-  4ceJf2s2  sin  y  4 


ccss2  sin  y  107ces28in  y  4cK6  sin  y 


36 


2cAVsin  y  23cessin  y  +  28ce  sin  y  4c/^9Cos  y 


360 


9 


5ce2s2 sin(20)  5es2sin(20)  5ee2cos(2u/)0  5ee2s2  cos(2w)0 

+  16  +  12  75s2  -  60  +  8 


ceJcos(2u>)0 


12 


4-  5cKjs  8  — 


5  ce2s20  5 es20  ce28  cO 


4 


6  +  2  ' 


(BA) 


24  3 

Substitution  of  (6.2)  -  (6.5),  (6.10),  (6.16),  (6.17),  (B.3),  and  (B.4)  into  equation 
(5.24)  results  in 


d2u2 

Ho2 


4  u2  =  2e  cos  y  y2  — 


120e4  cos(4y  —  20)  5 e4s4  cos(4y  —  20) 


4 


36000s2  -  28800 

7eK*,s 2  cos(4y  —  20)  31e4s2  cos(4y  -  20) 

.  _ 

e4  cos(4y  —  40) 


16 


4 


—  3e/f5  cos(4y  —  20) 
9e4  cos(4y  —  40) 


4 


4 


4 


e4cos(4y  — 20) 

240  +  15000s4  -  24000s2  4  9600  '  36000s2  -  28800 

29e4s4  cos(4y  -  40)  3 eK5s2  cos(4y  —  40)  923e4s2  cos(4y  -  40) 


192 


5760 


eK5  cos(4y  —  40)  89e4  cos(4y  —  40)  3e3s4  cos(3y  4  40) 


7200 


16 


e3s4  cos(3y  4  20)  21e3s2 cos(3y  4  20)  29e3  cos(3y  4  20) 

""48  ”  +  8  12 


480e3cos(3y  -  20)  5ess4  cos(3y  -  20) 
36000s2  -  28800  16 

5e3s2  cos(3y  -  20) 

— - - - -  4  8 Ks  cos(3y  -  20)  4 


—  10K$s2  cos(3y  -  20) 

l7e3cos(3y  -  20) 

30 


4 


35e2s4  cos(2y  4  40)  3eK7s2  cos(2y  4  20  4  2w) 

32  +  4 

3? Kss2  cos(2y  4  20  -  2w)  360e4  cos(2w)  cos(2y  4  20) 

4  +  36000s2  -  28800 

5e4s4  cos(2u>)  cos(2y  +  20)  19e4s2  cos(2w)  cos(2y  4  20) 


16 


96 


100 


+  e*  cos(2u;)  cos(2y  +  20)  _  5c2/f2s4cos(2y  +  20)  95 c2s4  cos(2y  +  26) 


+ 

+ 

+ 


80  2  12 
5 e2K2s2  cos(2y  +  26)  65 e2s2  cos(2y  +  26)  3tK7s 2  cos(2y  -26  +  2uj) 

4  8  +  4 

ZtKjS2  cos(2y  -  26  -  2u>)  4e4  cos(2u>)  cos(2y  —  20) 

4  15000s4  -  240005*  +  9600 

172c4  cos(2u>)  cos(2y  -  20)  e4s4  cos(2w)  cos(2y  -  20) 

36000s2  -  28800  8 


7eV  cos(2u>)  cos(2y  -  20)  23c4  cos(2u>)  cos(2y  -  20) 

80  +  3600 

576e2Kj  cos(2y  -  20)  164e4  cos(2y  -  20)  2400c2  cos(2y  -  20) 

36000s2  -  28800  36000s2  -  28800  +  36000s2  -  28800 

„2  is  e*s*  cos(2y  -  20)  257c2s4  cos(2y  -  20) 


—  e2K7sA  cos(2y  -  20)  - 


72 


5eK7s2  cos(2y  —  20)  13e*r5s2cos(2y  -  20)  I7e2/f2s2  cos(2y  -  20) 

2  +  6  +  30 

29e4s2  cos(2y  -  20)  245c2s2  cos(2y  -  20) 

+  ^5 - - 2eKTcos(2y  -  20) 

5eK$  cos(2y  -  20)  e2/f2  cos(2y  -  20)  67e4  cos(2y  -  20) 

3  +  50  +  3600 


+ 


+ 


e^cospy  -  20)  6c4  cos(2y  -  40)  528e2  cos(2y  -  40) 

12  36000s2  -  28800  36000s2  -  28800 

e4s4  cos(2y  -  40)  5e2s4  cos(2y  -  40)  3eK7s 2  cos(2y  -  40) 

2  16  ~  T 

9eK*,s 2  cos(2y  -  40)  121e4s2  cos(2y  -  40)  23e2s2  cos(2y  -  40) 


240 


240 


3cK~5  cos(2y  -  40)  29e4  cos(2y  -  40)  lie2  cos(2y  -  40) 


3 eK7s2  cos  (2 y  2u:) 


800 


+  eK7  cos(2y  +  2w)  - 


600 


3eK5s2  cos(2y  -  2 u>) 


eK5  cos(2y  -  2w)  - 


7e3s4 cos(y  +  40)  7es4  cos(y  +  40) 


8 


7K7s2  cos(y  ^  20  ^  2u)  7K5s2cos(y  ^  20  -  2u;) 


101 


+ 


+ 


2240c3  cos(2w)  cos(y  +  25)  5ess4  cos(2u>)  cos(y  +  26) 
36000s2  -  28800  4 

7ess2  cos('  w)  cos(y  +  26)  7escos(2i*>)  cos(y  +  26) 

12  +  90 


-  IOCS'S4  cos(y  +  26)  -f 


661css4  cos(y  +  26)  Alts*  cos(y  +  26) 


144 


+ 


36 


+ 

+ 


+ 


+ 


lOeKjs7  cos(y  +  26)  131ess1  cosfy  +  26)  ,  .. 

- 2 — - L - - L  _  5es.  cos(y  +  26) 

COSj'y  —  ^  +  2K7s2  cos(y  -2 6  +  2u)  +  2 Kisi  cosfy  -  26  -  2u>) 

3b 

480e3 cosfy  —  26)  25ess4  cosfy  —  26)  ,,  ,  , 

36000s2  -  28800  16  vy  ’ 

n  3  ,  ,  es2cos(y-25)  .  r .  .  17cs cosfy  -  26) 

2 eV  cos(y  -  25) - — - -  -  8tf7  cos(y  -  25)  + - — - 

6  30 

1120t3cos(y  —  45)  5ess4  cos(y  —  45)  53s4  cos(y  —  45) 

’36000s2  -  28800  +  48  12 

5 K7s2  cos(y  —  45)  7ess2  cos(y  -  45)  7es  cos(y  —  45)  570e4  cos(4y) 

3  24  180  36000s2  -  28800 

5e4s4cos(4y)  5e/sf5s2  cos(4y)  19e4s2  cos(4y)  5eK5cos(4y) 

32  +  4  +  192  2 

13e4cos(4y)  1120c3  cos(3y)  269ess4  cos(3y)  5K$s2  cos(3y) 

80  36000s2  -  28800  48  3 

679e3s2  cos(3y)  677escos(3y)  816e4  cos(2u>)  cos(2y) 

72  180  +  36000s2  -  28800 

lle4s4  cos(2w)  cos(2y)  359e4s2  cos(2w)  cos(2y) 

8  240 

17e4  cos(2a>)  cos(2y)  6e4cos(2y)  528c2cos(2y) 


600 


11  e2K2s*  cos(2y)  + 


36000s2  -  28800  36000s2  -  28800 

e4s4cos(2y)  695e2s4  cos(2y)  7 eK7s2  cos(2y) 

96 


2  96  4 

5eK5s2  cos(2y)  21e7K2s:  cos(2y)  121c4s2  cos(2y) 

4  2  240 

3007e2s2  cos(2y)  e/C5cos(2y)  29e4cos(2y) 

240  +  6  7C0S'  +  2  +  800 
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1137c2  cos(2y)  +  960es  cos(2u>)  cos  y  ^  15esa4cos(2u;)  cos  y 


200 


36000s2  -  28800 


15ess2  cos(2w)  cos  y  cscos(2u;)  cos  y  ,,  . 

-  - 7 — - -  + - — - - -  +  30eif2s4  cos  y 

4  oO 

15ess4  cos  y  275es4  cos  y  ^  ,,  ,  3ess2  cos  y 

+ - — - -  -  26eifjS2  cos  y - - 


16 


24 


61es2cos  y  7c8 cos  y  5eKjS7  cos(20  +  2w)  e-fos2  cos(20  +  2u>) 
_  - - - + - - - + - - - + - - - 

eKjs7  cos(20  —  2w)  5eJf5s2  cos(20  -  2w)  18Oe4cos(40) 

+  4  +  4  36000s2  -  28800 


+ 


5e4s4cos(40)  65e2s4  cos(40)  5s4cos(40)  5eK7s7  cos(40) 

8  4 


32 


32 


+ 


19e4s2 cos(40)  e4cos(40)  432e4  cos(2u>)  cos(20) 

192  160  +  36000s2  -  28800 

1056e2  cos(2w)  cos(20)  3c4s4  cos(2w)  cos(20) 

36000s2  -  28800  8 

e2s4  cos(2cl>)  cos(20)  19e4s2  cos(2u;)  cos(20) 

2  +  80 

23e2s2  cos(2w)  cos(20)  3e4cos(2w)  cos(20)  lle2cos(2u;)  cos(20) 
120  +  200  +  300 

276e4  cos(20)  „  ,  , 4  ln..  AT.  .  ...  7eVcos(20) 

36000,-  -  28800  “  *’*.**“•(">  "  <**‘‘“•<“1 - ^ 

127c2s4  cos(20)  7s4cos(20)  9eK7s2  cos(20)  ^  3e2K2s2  cos(20) 


16 


■+■  ,/C2s2  cos(20)  + 


6  2  2 
193e4s2 cos(20)  95e2s2 cos(20)  2s2cos(20) 


480 


T  _  ....  19e4cos(20) 

-  4eK7cos(20)  H - - -  -t- 


12 

2e4  cos(2w)2 


50e4  cos(2u;)2 


300  15000s4  -  24000s2  +  9600  36000s2  -  28800 

e4s4  cos(2w)2  23c4s2  cos(2cj)2  7c4cos(2w)2  576e2X2  cos(2w) 


32 


960 


3600 


36000s2  -  28800 


304e4  cos(2u.’)  2400e2  cos(2w)  c2/C2s4  cos(2w)  17e4s4  cos(2u-') 


36000s2  -  28800  36000s2  -  28800  '  2 

5e2s4  cos(2u;)  3eK7s2  cos(2w)  3eK5s2  cos(2w) 


24 
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l9eiKis2  cos(2u>) 
60 


533e4s2  cos(2cj) 
720 


41e2s2  cos(2u>) 
24 


+  cKjcos(2uj ) 


+ 

+ 

+ 

+ 


eKs  cos(2u;)  — 


e.2Ki  cos(2w)  19e4  cos(2u;) 


50 


+ 


1800 

17c4 


15000s4  -  24000s2  +  9600  36000s2  -  28800 

7eK 7s2 


25e4s4  437 e2s4  97s4 

+ 


192 

139eV 

72 


288 


8 


12 


-  K\s 2  - 


„  2  cK7  17c4 
”  1  S  +  3  +  720 


lie2 


+  2, 


c2  cos(2u;) 

12 

+  2 K\s*  +  17g  KzS  +  20 K2s4 
3 

31e2A'js2  2  847e4s2 

— “  UK‘S  -  “5760” 

(3.5) 


l 
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APPENDIX  C 


THE  CRITICAL  INCLINATION 


This  appendix  contains  a  discussion  of  the  critical  inclination  to  include  its 
mathematical  and  physical  basis  and  the  methods  used  by  various  authors  to  deal 
with  the  problem. 

In  Chapter  6  it  was  found  that  the  divisor  j  often  appeared  in  sec¬ 

ond  order  terms.  A  brief  explanation  was  given  for  the  appearance  of  this  divisor, 
and  it  was  shown  that  the  apparent  singularity  may  be  dealt  with  in  this  analy¬ 
sis  in  a  straight-forward  fashion.  However,  in  other  methods  the  problem  divisor 
has  caused  enormous  complications.  This  divisor  or  “critical  inclination”  problem 
appears  to  be  universal.  Hughes’  investigations  [Ref.  32]  reveal  that  every  ana¬ 
lytical  theory  on  orbit  perturbation  analysis  contains  the  critical  inclination  as  an 
inherent  characteristic.  As  stated  in  Chapter  2,  the  divisor  causes  an  irremovable 
singularity  in  Kozai’s  method.  Its  presence  has  prompted  many  authors  to  devise 
very  elaborate  techniques  to  get  around  the  problem  with  the  result  that  many 
otherwise  elegant  methods  are  somewhat  diminished. 

The  mathematical  reason  for  the  appearance  of  the  problem  term  is  as  fol¬ 
lows.  It  is  recalled  that  one  of  the  secular  effects  of  the  Earth’s  equatorial  bulge  is 
to  cause  the  line  of  apsides  to  rotate.  A  first  order  approximation  for  this  rotation 
is  given  by  the  rate  of  change  of  the  argument  of  perigee  which  from  Allan  [Ref. 


33  is 


cL-  3J2n  R‘  2. 


It  is  readily  seen  that  if  i  =  63°26',  the  perigee  position  is  stationary  and  the 


line  of  apsides  does  not  rotate.  When  the  usual  Poisson  method  of  successive 
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approximations  in  terms  of  a  small  parameter  is  adopted  for  perturbation  analysis, 
the  above  equation  appears  in  the  denominator  of  higher  order  terms. 

Coffey  [Ref.  34]  gives  an  excellent  background  discussion  on  the  history  of 
the  critical  inclination  problem.  He  states  that  A.  A.  Orlov  (1957)  was  perhaps 
the  first  to  notice  the  unusual  situation  at  the  critical  inclination.  Orlov  found 
that  the  conventional  Lindstedt- Poincare  algorithm  failed  at  this  critical  inclina¬ 
tion  and  at  zero  inclination.  Nevertheless,  Orlov  ignored  these  exceptional  cases 
and  developed  his  theories  without  addressing  the  problem  of  singularities.  Krause 
(1952),  Roberson  (1957),  and  Herget  and  Musen  (1958)  had  all  overlooked  this 
difficulty  in  their  assessment  of  the  long  term  effects  of  the  J2  term  on  the  Kep- 
lerian  ellipse  before  Brouwer  (1958)  finally  pointed  out  the  problem  to  the  latter 
authors.  Brouwer  hoped  that  his  use  of  von  Ziepel’s  method  of  eliminating  vari¬ 
ables  through  canonic  transformations  would  allow  him  to  dispose  of  the  terms 
that  lead  to  singularities.  However,  as  was  pointed  out  in  Chapter  2,  his  solution 
contained  singularities  at,  not  only  the  critical  inclination,  but  also  zero  inclination 
or  eccentricity. 

Finding  no  method  to  fit  the  critical  inclination  problem  into  current  theories, 
astronomers  have  devised  alternative  theories  to  deal  with  orbital  inclinations  near 
63°. 4.  The  critical  inclination,  they  reasoned,  was  a  small  divisor  problem  and 
therefore  could  be  handled  in  the  standard  way  introduced  by  Bohlen  [Ref.  35] . 
That  is,  in  the  vicinity  of  the  critical  inclination,  the  solution  may  be  expressed 
in  terms  of  the  square  root  of  a  small  parameter,  i.e.,  \/~J  instead  of  J.  [Ref.  36], 
The  resulting  solutions  are  totally  different  in  character  to  the  non-singular  case 
[Ref.  30].  For  example,  instead  of  precessing  secularly  the  argument  of  perigee 
will  tend  to  oscillate  about  a  central  position.  Taff  notes  the  lack  of  mathematical 
rigor  in  this  method  of  handling  the  critical  inclination  by  stating  that  expansions 
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in  y/1  and  theories  about  the  libration  of  perigee  are  a  result  of  misapplication  of 
perturbation  theories  [Ref.  18,  p  340j. 

Usually  when  singularities  arise  in  celestial  mechanics  there  is  a  physical 
reason  for  their  occurrence.  For  example,  singularities  in  the  solutions  for  the 
motion  of  asteroids  and  natural  satellites  are  caused  by  the  fact  that  their  orbital 
periods  are  nearly  commensurate  with  that  of  the  disturbing  body.  Under  the  same 
circumstances,  an  artificial  satellite  may  encounter  resonant  perturbations  in  its 
orbital  elements  caused  by  tesseral  harmonics.  However,  such  a  physical  reason 
does  not  seem  to  suggest  itself  for  the  critical  inclination  problem.  Message  [Ref. 
37]  implied  that  there  may  be  a  resonance  between  the  satellite’s  mean  motion 
relative  to  perigee  and  its  mean  motion  relative  to  the  ascending  node;  however, 
most  skeptics  have  not  been  convinced  [Ref.  32].  If  there  were  resonance  one  would 
expect  to  be  able  to  experimentally  measure  deviations  caused  by  resonances, 
but  satellite  guidance  engineers  have  reported  no  adverse  effects  at  the  critical 
inclination. 

Another  check  on  the  physical  effects  of  resonance  should  be  a  numerical 
check  of  the  equations  of  motion  at  the  critical  inclination;  however,  this  too 
has  been  inconclusive.  Lubowe  [Ref.  38]  carried  out  a  study  in  which  he  inte¬ 
grated  a  satellite’s  equations  of  motion  for  periods  up  to  24  hours  with  various 
initial  conditions  at  inclinations  near  and  away  from  the  critical.  He  found  no 
discernible  features  in  the  perturbations  which  distinguished  the  critical  from  any 
other.  However,  in  a  separate  analysis,  Hughes  came  to  the  opposite  conclusion. 
Hughes  developed  the  Hamiltonian  in  terms  of  the  Hill  variables  and  showed  that 
the  critical  inclination  remained  an  inherent  part  of  the  solution  throughout  a 
multitude  of  various  canonic  transformations.  He  then  integrated  the  equations 
of  motion  numerically  and  was  able  to  show  some  resonant  effects  in  the  perigee 
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height  for  satellites  at  the  critical  inclination.  He  then  reasoned  from  this  analysis 
that  the  critical  inclination  is  unique,  that  it  represents  actual  physical  resonance 
and  that  it  is  not  merely  a  by-product  of  the  method  of  solution  or  the  type  of 
variable  used  in  the  analysis. 

This  disagreement  characterizes  the  lingering  debate  over  the  question  of 
whether  the  critical  inclination  is  an  artifact  of  the  analysis  process  or  a  physical 
reality.  Most  textbooks  choose  to  mention  the  problem;  however,  few  voice  an 
opinion.  Roy  [Ref.  20]  devotes  only  two  sentences  to  the  subject,  merely  noting 
that  some  perturbation  analyses  break  down  there.  Hagihara  [Ref.  39],  mentions 
that  the  critical  inclination  is  an  essential  feature  of  the  perturbation  process 
but  notes  that  heated  discussions  continue  concerning  its  physical  reality.  Other 
authors,  namely  Herrick  [Ref.  40],  Kovalevsky  [Ref.  41],  and  Geyling  [Ref.  li 
mention  the  problem  and  the  standard  procedure  for  dealing  with  it,  but  they 
make  no  mention  of  physical  effects.  Taff  [Ref.  18,  p  340]  takes  a  firmer  position 
by  stating  that  it  is  not  a  prediction  of  second-order  perturbation  theory  that 
there  are  infinitely  large  or  infinitely  rapid  changes  in  the  orbital  elements,  and 
that  the  correct  resolution  of  these  “unphysicar  predictions  is  not  to  rely  on  bad 
mathematics. 

This  present  analysis  lends  support  to  Taff’s  assertion.  Although  the  criti¬ 
cal  inclination  problem  was  manifested  in  the  perturbation  process,  it  was  shown 
that  it  does  not  cause  singularities,  nor  are  there  any  unusual  physical  effects  pre¬ 
dicted  by  this  theory.  In  fact,  it  was  a  remarkable  aspect  of  this  analysis  that 
in  the  limit  as  the  inclination  approached  the  critical,  all  potential  singularities 
were  canceled.  This  is  not  true  of  the  analyses  which  use  the  technique  of  De- 
launey  /von  Ziepel  canonic  transformations.  A  remarkable  aspect  of  the  canonic 
transformation  method  is  that  the  non-integrability  of  the  equations  of  motion 


is  not  obvious  when  the  system  is  expressed  in  coordinates  or  elements,  but  after 
the  canonical  transformation  it  is  quite  clear  [Ref.  42].  The  critical  inclination  is 
an  intrinsic  singularity  of  the  method  [Ref.  32]. 
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